Dispersion extraction for acoustic data using time frequency analysis

ABSTRACT

This invention pertains to the extraction of the slowness dispersion characteristics of acoustic waves received by an array of two or more sensors by the application of a continuous wavelet transform on the received array waveforms (data). This produces a time-frequency map of the data for each sensor that facilitates the separation of the propagating components thereon. Two different methods are described to achieve the dispersion extraction by exploiting the time frequency localization of the propagating mode and the continuity of the dispersion curve as a function of frequency. The first method uses some features on the modulus map such as the peak to determine the time locus of the energy of each mode as a function of frequency. The second method uses a new modified Radon transform applied to the coefficients of the time frequency representation of the waveform traces received by the aforementioned sensors. Both methods are appropriate for automated extraction of the dispersion estimates from the data without the need for expert user input or supervision.

FIELD OF THE INVENTION

This invention is generally related to acoustic data analysis, and moreparticularly to processing acoustic waveforms where there is dispersion(frequency dependent wave velocity and attenuation), including but notlimited to acoustic analysis of materials and subsurface formations.

BACKGROUND OF THE INVENTION

Dispersive processing of borehole acoustic data is useful for thecharacterization and estimation of rock properties using boreholeacoustic data containing propagating modes. It is both more complex andmore informative than non-dispersive processing, where bulk velocitiesin the rock are measured without any frequency dependence. Perhaps themost common parameters used to describe the dispersion characteristicsare the wavenumber, k(f), and the attenuation, A(f), both of which arefunctions of the frequency f and are of interest in the oil industry forcharacterizing the properties of reservoirs. The dispersioncharacteristics include the group and phase slowness (reciprocal ofvelocity) as a function of frequency linked to the wavenumber, k(f), asfollows:

$\begin{matrix}{{{s_{\phi}(f)} = {\frac{1}{V_{phase}(f)} = \frac{k(f)}{f}}}{and}} & (1) \\{{s_{g}(f)} = {\frac{1}{V_{group}(f)} = \frac{\mathbb{d}{k(f)}}{\mathbb{d}f}}} & (2)\end{matrix}$These two quantities are not independent, and satisfy:

$\begin{matrix}{s_{g} = {s_{\phi} + {f{\frac{\mathbb{d}s_{\phi}}{\mathbb{d}f}.}}}} & (3)\end{matrix}$

One known class of dispersive processing solutions uses physical modelswhich relate the rock properties around the borehole to predicteddispersion curves. Waveform data collected by an array of sensors isbackpropagated according to the modeled dispersion curves, and the modelis adjusted until there is good semblance among these backpropagatedwaveforms, thereby indicating a good fit of the model to the data.However, models are presently available for only simpler cases, andother physical parameters still need to be known. Further, there is arisk that biased results will be produced if there is model mismatch orinput parameter error. Moreover, this class of solutions assumes thepresence of only a single modeled propagating mode, and pre-processingsteps such as time windowing and filtering may be desireable to isolatethe mode of interest. These also require user input, and in some casesthe results may be sensitive to the latter, requiring expert users forcorrectly processing the data.

One way to mitigate these drawbacks is to directly estimate thedispersion characteristics from logging data, i.e., without reference toparticular physical models. Not only can this be used for quantitativeinversion of parameters of interest, but the dispersion curves carryimportant information about the acoustic state of the rock and areimportant tools for interpretation and validation. Moreover, in order tobe part of a commercial processing chain, the dispersion estimationmethod should be capable of operation in an automated unsupervisedmanner and in particular the accuracy of the results should not dependon input from highly skilled users. One technique of estimatingdispersion characteristics directly from data collected by an array ofsensors, for example, in seismic applications, is to use a 2D FFT, alsocalled the f-k (frequency wavenumber) transform. This techniqueindicates the dispersion characteristics of propagating waves, bothdispersive and non-dispersive. However, its effectiveness is limited tolarge arrays of tens of sensors. For applications with fewer sensors,such as commonly found in sensor arrays (2-13 sensors) on boreholewireline and LWD tools, this technique may lack the necessary resolutionand accuracy to produce useful answers.

A high resolution method appropriate for shorter arrays was developedusing narrow band array processing techniques applied to frequencydomain data obtained by performing an FFT on the array waveform data.However, while this is an effective tool for studying dispersionbehavior in a supervised setting with user input, it produces unlabelleddots on the f-k plane frequency by frequency and therefore may not beentirely suitable for deployment as an automated unsupervised processingmethod. Moreover, the operation of the FFT washes out the informationpertaining to the time of arrival of various propagating modes, therebycompromising performance and the effectiveness of interpretation,especially for weaker modes of interest overlapping with stronger onesin frequency domain.

A parametric method for estimating dipole flexural dispersion is knownwhich can be used for automated dispersion extraction. However this islimited to the particular (flexural) dispersive mode, and cannot bereadily extended to the general case. Moreover, it is a one-componentextraction technique.

Recently, a new algorithm has been proposed using time frequencyanalysis with continuous wavelet transform along with beamformingmethods. However, the algorithm does not account for attenuation and hasa propagation model with a pure phase response. In real worldapplications, many modes of interest are attenuative and the attenuationis therefore a key parameter of interest. Further, the phase response isestimated using a 2-D search over delays and phase corrections appliedto the original data, for each of which, a time-frequency map isgenerated. Extending this approach to include attenuation would make ita 3-D search, which would be relatively computationally intensive forcommercial real time applications. Finally, the use of this approach aswell as the criterion of maximum energy for identifying propagatingmodes on the time-frequency maps generated above may obscure weak modesof interest by stronger in-band interfering modes even when they areseparated in time.

SUMMARY OF THE INVENTION

In accordance with one embodiment of the invention, a method ofdispersion extraction for acoustic data comprises the steps of:receiving a time series of acoustic data associated with a plurality ofsensors having different transmitter-receiver spacings; generating aplurality of time-frequency representations from the received timeseries; identifying a characteristic feature in said time-frequencyrepresentations; selecting at least one frequency; defining a kernelcomprising: estimating a time location of the characteristic feature atthe selected frequency; associating a corresponding characteristicfeature for ones of the spacings at the selected frequency; fitting atleast a segment of a line through at least some of the time locations ofthe associated characteristic features, plotted against thecorresponding spacings; estimating group slowness at the selectedfrequency based at least in-part on an indication of fitted line segmentslope; estimating a phase difference between sensors; estimating anattenuation between sensors; using the phase difference and theestimated group slowness to estimate phase slowness; and repeating thekernel for other frequencies of interest, obtaining a dispersion curvecomprising the group slowness, the phase slowness, and the attenuationas a function of frequency, and providing the dispersion curve intangible form.

In accordance with another embodiment of the invention, a method ofdispersion extraction for acoustic data comprises the steps of:receiving a time series associated with a plurality of sensors;generating a time-frequency representation for each sensor; selecting areference sensor; identifying a characteristic feature in the timefrequency representation for the reference sensor; estimating timelocations across frequency of the characteristic feature; selecting aninitial frequency at which energy is greater than a threshold value oversaid time locations; defining a kernel comprising: collecting timefrequency coefficients at the selected frequency for the sensors into arepresentation of data indexed by time and sensor; selecting timewindows centered on a set of positions in a local neighborhood of thetime location at the selected frequency; selecting a set of testmoveouts for ones of the time windows corresponding to expected range ofgroup slowness; shifting and aligning the representation of coefficientdata in the window corresponding to the test moveout and windowpositions; computing estimates of phase difference and attenuationacross the representation of shifted and aligned data for the testmoveouts and window positions; using the estimated phase difference andattenuation to perform a modified stacking operation on therepresentation of shifted and aligned data for the test moveouts andwindow positions to produce an output representation for a plurality ofvalues of time position and moveout; finding a maximum peak of theoutput representation and using the corresponding moveout to estimatethe group slowness; using a corresponding value of the computed estimateof the phase difference along with the group slowness to generate anestimate of the phase slowness; setting a corresponding value of thecomputed estimate of attenuation as an attenuation estimate; iterativelyincreasing frequency from the selected initial frequency comprising:selecting a higher frequency; constraining the moveout based on theprevious computed value of the group slowness; repeating the kernel forthe selected higher frequency until a final highest frequency is reachedor when the maximum computed semblance falls below a specifiedthreshold; iteratively decreasing frequency from the selected initialfrequency comprising selecting a lower frequency; constraining themoveout based on the previous computed value of the group slowness;repeating the kernel for the selected lower frequency until a finallowest frequency is reached or when the maximum computed semblance fallsbelow a specified threshold; obtaining a dispersion curve comprising thegroup slowness, the phase slowness, and the attenuation as a function offrequency; and providing the dispersion curve in tangible form.

In accordance with another embodiment of the invention, apparatus fordispersion extraction for acoustic data comprises: a sonic logging tooloperable to generate a time series of acoustic data associated with aplurality of sensors having different transmitter-receiver spacings; ananalyzer operable to: generate a plurality of time-frequencyrepresentations from the time series; identify a characteristic featurein said time-frequency representations; program code operable to selectat least one frequency; define a kernel, including: estimating a timelocation of the characteristic feature at the selected frequency;associating the corresponding characteristic feature for ones of thespacings at the selected frequency; fitting at least a segment of a linethrough at least some of the time locations of the associatedcharacteristic features;

estimating group slowness at the selected frequency based at leastin-part on slope of the fitted line segment; estimating a phasedifference between sensors; estimating attenuation between the sensors;using the phase difference and the group slowness to estimate the phaseslowness; and repeating the kernel for other frequencies of interest,obtaining a dispersion curve comprising the group slowness, the phaseslowness, and the attenuation as a function of frequency, and storingthe dispersion curve in memory.

In accordance with another embodiment of the invention, apparatus fordispersion extraction for acoustic data comprises: a sonic logging tooloperable to generate a time series associated with a plurality ofsensors; an analyzer operable to: generate a time-frequencyrepresentation for ones of the sensors; select a reference sensor;identify a characteristic feature in the time frequency representationfor the reference sensor; determine time locations across frequency ofthe said characteristic feature; select an initial frequency at whichenergy is maximum over the said time locations; define a kernel,including: collecting time frequency coefficients at the selectedfrequency for ones of the sensors into a representation of data indexedby time and sensor; selecting time windows centered on a set ofpositions in a local neighborhood of the time location at the selectedfrequency; selecting a set of test moveouts for ones of the time windowscorresponding to expected range of group slowness; shifting and aligningthe representation of coefficient data in the window corresponding toones of the test moveout and window positions; computing estimates ofphase difference and attenuation across the representation of shiftedand aligned data for ones of the test moveouts and window positions;using the estimated phase difference and attenuation to perform amodified stacking operation on the representation of shifted and aligneddata for ones of the test moveouts and window positions to produce anoutput representation for values of time position and moveout; finding amaximum peak of the output representation and using a correspondingmoveout to estimate group slowness; using a corresponding value of thecomputed estimate of the phase along with group slowness to generate anestimate of phase slowness; setting a corresponding value of thecomputed estimate of attenuation as an attenuation estimate; iterativelyincreasing frequency from the selected initial frequency comprisingselecting a higher frequency; constraining the moveout based on theprevious computed value of the group slowness; repeating the kernel forthe selected higher frequency until a final highest frequency is reachedor when the maximum computed semblance falls below a specifiedthreshold; iteratively decreasing frequency from the selected initialfrequency comprising selecting a lower frequency; constraining themoveout based on the previous computed value of the group slowness;repeating the kernel for the selected lower frequency until a finallowest frequency is reached or when the maximum computed semblance fallsbelow a specified threshold; obtaining a dispersion curve comprising thegroup slowness, the phase slowness, and the attenuation as a function offrequency; and providing the dispersion curve in tangible form.

The inventive technique advantageously allows the extraction of thegroup slowness dispersion along with the phase slowness dispersion andattenuation from the array waveform data without the use of any physicalmodel by exploiting the time frequency localization of the propagatingmodes at each sensor as well as the continuity of the dispersion curveas a function of frequency. In contrast with at least some prior arttechniques, in at least one embodiment of this invention thetime-frequency mapping is computed once and, using the continuity of thedispersion curves, a restricted 1-D search is utilized over delays overthe wavelet coefficients followed by closed form estimation of the bestcorresponding phase correction and attenuation. Moreover, the criterionof semblance that is effective in identifying weak modes so long as theyare non-overlapping in the time-frequency plane can be utilized. As aresult, an improved and practical means can be provided for directlyextracting dispersion characteristics from data collected by arelatively short array of sensors.

The invention has practical application for the processing of boreholeacoustic data dominated by dispersive borehole modes usinglogging-while-drilling and wireline tools, but could also be used insituations where dispersion exists and needs to be estimated. Thesecould include, for example, the processing of surface waves such asground roll in seismic data and the processing of dispersive Lamb modeswith ultrasonic data for non-destructive evaluation. These dispersioncharacteristics can then be used to extract information of interest suchas rock compressional and shear velocities, radial variation, anisotropyetc. In the Detailed Description which follows a method to estimate thedispersion characteristics (group and phase velocity and attenuation asa function of frequency) of one or more propagating modes using an arrayof two or more sensors will be described. While the examples providedherein are directed towards borehole applications, as described above,the methods and apparatuses described herein may also be applied to thecharacterization of other materials such as cement, and human or animaltissue.

Further features and advantages of the invention will become morereadily apparent from the following detailed description when taken inconjunction with the accompanying Drawing.

BRIEF DESCRIPTION OF THE DRAWING

FIG. 1 illustrates a sonic logging tool with an analyzer unit.

FIGS. 2 a and 2 b illustrate the efficacy of local Taylor expansions oforders one and two in capturing the behavior of a typical quadrupoledispersion curve in the wavenumber and slowness domains respectively.

FIG. 3 illustrates an example of propagation in a plane.

FIG. 4 illustrates a set of waveforms collected by an array of sensorsand the transformation of one of the sensor waveforms into a 2-D timefrequency map.

FIG. 5 illustrates the collection of 2-D maps obtained by thetransformation from time domain to time-frequency domain of eachwaveform from the sensor array.

FIG. 6 illustrates the computation of group slowness at one scale (i.e.,frequency).

FIG. 7 illustrates the result of using data association to track andlabel modulus peaks for various modes on the time frequency plane forthe reference sensor data.

FIG. 8 illustrates comparison of the dispersion extraction method(method 1) on a frame of oilfield data recorded in a fluid filledborehole, where the circles correspond to the dispersion curve extractedfrom the Prony's method while the cross marks correspond to thedispersion curve computed from the algorithm presented in this study.

FIG. 9 illustrates operation of the Exponential Projected RadonTransform (EPRT) on the complex coefficients of the CWT of array data ata particular scale a.

FIG. 10 illustrates the result of correcting for the bias in the groupslowness estimate for synthetic data.

FIG. 11 illustrates the use of the EPRT method to generate estimates ofphase slowness dispersion and attenuation on a frame of real oilfielddata compared with the results from the Prony method.

FIG. 12 illustrates a flowchart of the EPRT method.

DETAILED DESCRIPTION

In accordance with the invention, dispersion curves are produced by ananalyzer unit from acoustic data gathered by a borehole logging tool.Regardless of the embodiments illustrated, the logging tool may be of awireline type, logging while drilling type, or any other type. Furtherthe transmitter may be used to excite a dipole, monopole, quadrupole orany other multipole borehole acoustic mode. It should also be noted thatseismic and borehole seismic data may be utilized, since both cancontain dispersive arrivals. Each “sensor” includes a source (one ormore transmitters) and a receiver (one or more associated hydrophones).For example, one transmitter and eight receivers constitute a set ofeight sensors. Characteristic features of a sensor include the type ofacoustic disturbance emitted by the source, and the distance (spacing)between the source and receiver.

FIG. 1 illustrates one example of a logging tool (106) used to acquiredata from which dispersion curves will be extracted. The tool has aplurality of receivers and transmitters. The illustrated logging tool(106) also includes multi-pole transmitters such as crossed dipoletransmitters (120, 122) (only one end of dipole (120) is visible inFIG. 1) and monopole transmitters (109) (close) and (124) (far) capableof exciting compressional, shear, Stoneley, and flexural waves. Aplurality of receivers is arranged on the logging tool (106), spacedfrom the transmitters. In particular, the logging tool includes thirteenspaced receiver stations, each receiver station comprising multiplereceiver hydrophones (126) mounted azimuthally at regular intervalsaround the circumference of the tool. With thirteen receivers (126) ateach receiver station, there are a total of one-hundred and fourreceiver elements. Other configurations, such as a Digital Sonic Imaging(DSI) tool with four receivers at each of eight receiver stations, orincorporating other multipole sources such as quadrupole, are alsopossible. The use of a plurality of receivers and transmitters resultsin improved signal quality and adequate extraction of the variousborehole signals over a wide frequency band It is noted that thedistances and number of receivers and transmitters shown in this exampleare provided as exemplary only and are not intended to be limiting;other configurations may be employed.

During operation of the logging system (100), the subsurface formation(102) is traversed by the borehole (104) which may be filled withdrilling fluid or mud. The logging tool (106) is suspended from anarmored cable (108) and may have optional centralizers (not shown). Thecable (108) extends from the borehole (104) over a sheave wheel (110) ona derrick (112) to a winch forming part of surface equipment, which mayinclude an analyzer unit (114). Well known depth gauging equipment (notshown) may be provided to measure cable displacement over the sheavewheel (110). The tool (106) may include any of many well known devicesto produce a signal indicating tool orientation. Processing andinterface circuitry within the tool (106) amplifies, samples anddigitizes the tool's information signals for transmission andcommunicates them to the analyzer unit (114) via the cable (108).Electrical power and control signals for coordinating operation of thetool (106) may be generated by the analyzer unit (114) or some otherdevice, and communicated via the cable (108) to circuitry providedwithin the tool (106). The surface equipment includes a processorsubsystem (116) (which may include a microprocessor, memory, clock andtiming, and input/output functions—not separately shown), standardperipheral equipment (not separately shown), and a recorder (118). Thetool (106) also includes a set of instructions that, when executed,provide a series of dispersion curves extracted from data at each depth.The logging tool (106) is representative of any logging device that maybe used in accordance with principles described herein. It will beunderstood by those of skill in the art having the benefit of thisdisclosure that other suitable logging device, including LWD devices,can also be utilized.

Review of Complex Continuous Wavelet Transform

The continuous wavelet transform is a transformation allowing thedecomposition of an arbitrary time or space dependent signal, s(p), intoelementary contributions of functions called wavelets obtained bydilation and translation of a “mother” or analyzing wavelet g(p). Inthis document the terms “waveforms,” “signal,” “function,” and “timeseries” are used to refer to data collected by any of a set of receiversat a plurality of sampling points in time or space. Note that the datacan be viewed as a series (e.g., “time series”) that represents theevolution of the observed quantity as a function of time (or space),when plotted out versus time (or space), as tracing out the shape of theacoustic waves received (e.g., “waveform”), and also as containinginformation to be extracted (e.g., “signal”). For the purposes of thisdescription, let s(p) be an arbitrary time or space dependent signal andg(p) the chosen complex and progressive analyzing wavelet, necessary tostudy wave propagation phenomena, and let p be the time variable. Thecontinuous wavelet transform S(b,a) of a function s(p) is the scalarproduct of this signal by the dilated (contracted) and translatedwavelets family, g, such as

${T^{b}{D^{a}\left\lbrack {g(p)} \right\rbrack}} = {a^{{- 1}/2}{{g\left( \frac{p - b}{a} \right)}.}}$In the case of L₂ normalization this yields:

$\begin{matrix}{{{S\left( {b,a} \right)} = {\text{<}T^{b}{D^{a}\left\lbrack {g(p)} \right\rbrack}}},{{s(p)}>={a^{{- 1}/2}{\int{{s(p)}{\overset{\_}{g}\left( \frac{p - b}{a} \right)}{\mathbb{d}p}}}}}} & (4)\end{matrix}$while for L₁ normalization the expression can be represented as:

$\begin{matrix}{{{S\left( {b,a} \right)} = {\text{<}T^{b}{D^{a}\left\lbrack {g(p)} \right\rbrack}}},{{s(p)}>={a^{- 1}{\int{{s(p)}{\overset{\_}{g}\left( \frac{p - b}{a} \right)}{\mathbb{d}p}}}}}} & (5)\end{matrix}$where g is the complex conjugate of g dilated in time by a (a>0) andtranslated in time by b, homogeneous to the time in this case, (b ∈ R).a and b are respectively the scale (or dilation) parameter, that can beinterpreted like a zoom, and the translation parameter. Small dilationswill be related to the high frequencies and vice versa. In order tocorrectly define and give a physical meaning to the phase of the waveletcoefficients, the analyzing wavelet should satisfy the analytic orprogressive property (i.e.: ĝ(ω)=0, for negative (spatial or time)frequency components ω<0). The calculation of the wave fronts ofdifferent wave contributions and their spectral components can beperformed precisely without artifacts or interferences due to theabsence of Fourier components on the negative axis.

There exists some flexibility in the choice of the analyzing wavelet,but it should preferably satisfy the admissibility condition deducedfrom the isometric property of the transform in the following sense:there exists for every s(t) a constant c_(g) depending only on thewavelet g such that:

$\begin{matrix}{{{\int{{{s(t)}}^{2}{\mathbb{d}t}}} = {c_{g}^{- 1}{\int{\int{{{S\left( {b,a} \right)}}^{2}\frac{{\mathbb{d}a}{\mathbb{d}b}}{a^{2}}}}}}},{and}} & (6) \\{c_{g} = {{2\;\pi\;{\int{\frac{{{\hat{g}(\omega)}}^{2}}{\omega }{\mathbb{d}\omega}}}} < {\infty.}}} & (7)\end{matrix}$ĝ is the Fourier transform of g with ω as the dual variable of the timet and the inequality on the right is the admissibility condition. Itfollows that g is of zero mean (∫g(t)dt=0 or ĝ(0)=0 ). If this conditionis satisfied, there exists an inversion formula which reconstructs theanalyzed signal (Grossmann, A. and Morlet, J., 1984, Decomposition ofHardy functions into square integrable wavelets of constant shape,SIAM—J. Math. Anal., 15, 723-736).

$\begin{matrix}{{s(t)} = {{Re}\left\lbrack {c_{g}^{- 1}{\int{\int{{S\left( {b,a} \right)}a^{{- 1}/2}{g\left( \frac{t - b}{a} \right)}\frac{{\mathbb{d}a}{\mathbb{d}b}}{a^{2}}}}}} \right\rbrack}} & (8)\end{matrix}$where Re[.] represents the real part.

It can be shown that s(t) can also be obtained by a simple inversionformula involving only a one-dimensional integral over the dilationparameter:

$\begin{matrix}{{{s(t)} = {{Re}\left\lbrack {k_{g}^{- 1}{\int{{S\left( {t,a} \right)}\frac{\mathbb{d}a}{a^{3/2}}}}} \right\rbrack}},} & (9)\end{matrix}$provided that the analyzing wavelet satisfies the followingadmissibility condition:

$\begin{matrix}{{{\int{\frac{{\hat{g}(\omega)}}{\omega }{\mathbb{d}\omega}}} < \infty},{and}} & (10) \\{k_{g} = {2\;\pi{\int{\frac{\hat{g}(\omega)}{\omega}{{\mathbb{d}\omega}.}}}}} & (11)\end{matrix}$Since the CWT is non-orthogonal, <g(b,a), g(b′,a′)>≠0. There exists areproducing kernel N_(g) defined from equations (4) and (7) as:N _(g)(b,a,v,u)=c _(g) ⁻¹ <g(b,a),g(v,u)>.   ((12)LetN _(g)(v,u;0,1)=c _(g) ⁻¹ <T ^(v) D ^(u) g,g>.   (13)N_(g) verifies:

$\begin{matrix}{{N_{g}\left( {v,{u;b},a} \right)} = {{N_{g}\left( {0,{1;\frac{b - v}{u}},\frac{a}{u}} \right)}.}} & (14)\end{matrix}$N_(g) is then a function of scale ratio

$\frac{a}{u}$with a distance defined in scale ratio

$\frac{b - v}{u}.$It reflects the coherence of the wavelet coefficient in the half-plane(b,a) (time-scale resolution of the wavelet). Its modulus |N_(g)| has amaximal value when the pairs {b,a}={v,u} and decays quickly when thedistance of {b,a} from {v,u} increases, it means

${\frac{a}{u} \approx 1},{\frac{b - v}{u}{{\bullet 1}.}}$

A progressive analyzing wavelet such as a Morlet type in which

${{g(t)} = {{\exp\left( {{\mathbb{i}\omega}_{0}t} \right)}{\exp\left( \frac{- t^{2}}{2\beta^{2}} \right)}}},$with ω₀=2π, β=1 yielding ĝ(ω)≈0 for ω<0 is selected. It is recognizedthat the Morlet wavelet is not a true wavelet in that its integral isnot zero. However, for a large enough ω₀ (in practice larger than 5),the integral of the Morlet wavelet is small enough that it can be usednumerically as if it were a wavelet (Grossmann, A., Kronland-Martinet,R., Morlet, J., 1989, Reading and understanding continuous wavelettransform. Wavelet, Time-frequency Methods and Phase Space, Ed. J. M.Combes, A. Grossmann, P. Tchamitchian, Springer-Verlag, Berlin).Denoting

${a^{\prime} = {{\frac{a}{u}\mspace{14mu}{and}\mspace{14mu} b^{\prime}} = \frac{b - v}{u}}},$and using results from Gradshteyn, I. S. and Ryzhik, I. M., 1990, Tableof Integrals, Series and Products, Academic Press, New York, the modulusand the phase of the reproducing kernel has the explicit form:

$\begin{matrix}{{{{N_{g}\left( {0,{1;b^{\prime}},a^{\prime}} \right)}} = {\frac{\beta}{c_{g}}\sqrt{\frac{2\;\pi\; a^{\prime}}{1 + a^{\prime 2}}}{\exp\left( {{- \frac{1}{2}}\frac{\frac{\begin{matrix}{{\omega_{0}{\beta^{2}\left( {a^{\prime} - 1} \right)}^{2}} +} \\b^{\prime 2}\end{matrix}}{\beta^{2}}}{1 + a^{\prime 2}}} \right)}}}{{\arg\left\lbrack {N_{g}\left( {0,{1;b^{\prime}},a^{\prime}} \right)} \right\rbrack} = \frac{\omega_{0}{b^{\prime}\left( {a^{\prime} + 1} \right)}}{1 + a^{\prime 2}}}} & (15)\end{matrix}$It follows that:

$\begin{matrix}{{{{{N_{g}\left( {v,{u;b},a} \right)}} = {\frac{\beta}{c_{g}}\sqrt{\frac{2\;\pi\;{ua}}{u^{2} + a^{2}}}{\exp\left( {{- \frac{1}{2}}\frac{\begin{matrix}{{\omega_{0}{\beta^{2}\left( {a - u} \right)}^{2}} +} \\\frac{\left( {v - b} \right)^{2}}{\beta^{2}}\end{matrix}}{u^{2} + a^{2}}} \right)}}};}{{\arg\left\lbrack {N_{g}\left( {v,{u;b},a} \right)} \right\rbrack} = \frac{{\omega_{0}\left( {b - v} \right)}\left( {u + a} \right)}{u^{2} + a^{2}}}} & (16)\end{matrix}$

From equations (4) and (7), wavelet coefficients satisfy the reproducingequation:

$\begin{matrix}{{S\left( {v,u} \right)} = {\int{{S\left( {b,a} \right)}{N_{g}\left( {v,u,b,a} \right)}{\frac{{\mathbb{d}b}\;{\mathbb{d}a}}{a^{2}}.}}}} & (17)\end{matrix}$This allows the use of the interpolation formula introduced by A.Grossmann (Grossmann, A., Kronland-Martinet, R., Morlet, J., 1989,Reading and understanding continuous wavelet transform. Wavelet,Time-frequency Methods and Phase Space, Ed. J. M. Combes, A. Grossmann,P. Tchamitchian, Springer-Verlag, Berlin), to reconstruct an approximatevalue of the CWT from the value of the Discrete Wavelet Transform (DWT).

Wavelet Transform Local Modulus Maxima

As will be described in greater detail below, the local maximum (pluralmaxima) of the wavelet transform modulus can be used to implement theinventive technique. For the purposes of this description, local extremaof the wavelet transform of s(t), is any point (a₀,t_(o)) such that

$\begin{matrix}{\frac{\partial{S\left( {a_{0},t_{o}} \right)}}{\partial t} = 0.} & (18)\end{matrix}$The local maxima of the wavelet transform of s(t), is any point(a₀,t_(o)) such that when t belongs to either the right of the leftneighborhood of t_(o), |S(a₀,t)|<|S(a₀,t₀)| and when t belongs to theother side of the neighborhood of t_(o), |S(a₀,t)|≦|S(a₀,t₀)|. Themaxima line of the wavelet transform is any connected curve in the scalespace (a,t) along which the points are local maxima of the wavelettransform. A local maximum (a₀,t_(o)) of the wavelet transform isstrictly maximum either on the right or the left side of the t₀. Itshould be noted that the term “local maximum” of the wavelet transformmay be used to refer to a local maximum of the wavelet transform modulusin order to simplify the explanation.

Reassigned Scalogram

A “reassigned scalogram” approach may be useful to separate componentswhich are too close in the time-scale plane by pre-processing, i.e.,before applying the dispersion extraction methods described below. Abrief description is given here, and a more general description can befound in Auger, F., and Flandrin, P., 1995, Improving the readability oftime frequency time-scale representation by the reassignment method:IEEE Trans. Sig. Proc., 43, No. 5, 1068-1089. The reassignment principlecan be applied to time-scale representations of the affine class (O.Rioul and P. Flandrin, Time-scale energy distributions: a general classextending wavelet transforms, IEEE Trans. Signal Process., SP-40(7),1746-1757, 1992). A widely used member of this class is the scalogram,which is the squared modulus of the continuous wavelet transform:

$\begin{matrix}{{{{SC}_{s}^{g}\left( {t,a} \right)} = {{{W_{s}^{g}\left( {t,a} \right)}}^{2}\mspace{14mu}{with}}}{{W_{s}^{g}\left( {t,a} \right)} = {\frac{1}{\sqrt{a}}{\int_{- \infty}^{+ \infty}{{s(u)}{\overset{\sim}{g}\left( \frac{u - t}{a} \right)}\ {\mathbb{d}u}}}}}} & (19)\end{matrix}$where g(t) is the mother wavelet described previously. The scalogramresults from an affine smoothing of the Wigner Ville distribution and isdefined as:

$\begin{matrix}{{{SC}_{s}^{g}\left( {t,a} \right)}{\int{\int{{W_{s}\left( {v,\xi} \right)}{W_{g}\left( {\frac{v - t}{a},{a\;\xi}} \right)}\frac{{\mathbb{d}v}{\mathbb{d}\xi}}{2\;\pi}}}}} & (20)\end{matrix}$where W_(s) (□,□) and W_(g) (□,□) denote the Wigner-Ville distributionsof s and g respectively. As evidenced by this expression, SC_(s) ^(g)(t,a) can be interpreted as the summation of a whole set of energy measuresW_(s)(v, ξ) contained within a time frequency domain delimited by

${W_{g}\left( {\frac{v - t}{a},{a\;\xi}} \right)}.$Instead of assigning this number to the geometric center of this domain,which does not depend on the analyzed signal, it may be assigned to thecenter of gravity, defined by:

$\begin{matrix}\begin{matrix}{{\hat{t}\left( {t,a} \right)} = {\frac{1}{{SC}_{s}^{g}\left( {t,a} \right)}{\int{\int{{{vW}_{s}\left( {v,\xi} \right)}{W_{g}\left( {\frac{v - t}{a},{a\;\xi}} \right)}\frac{{\mathbb{d}v}{\mathbb{d}\xi}}{2\;\pi}}}}}} \\{{\hat{\omega}\left( {t,a} \right)} = \frac{\omega_{o}}{\hat{a}\left( {t,a} \right)}} \\{= {\frac{1}{{SC}_{s}^{g}\left( {t,a} \right)}{\int{\int{\xi\;{W_{s}\left( {v,\xi} \right)}{W_{g}\left( {\frac{v - t}{a},{a\;\xi}} \right)}\frac{{\mathbb{d}v}{\mathbb{d}\xi}}{2\;\pi}}}}}}\end{matrix} & (21)\end{matrix}$The resulting reassigned scalogram, defined as:

$\begin{matrix}{{\overset{\Cup}{S}{\overset{\Cup}{C}}_{s}^{g}} = {\int{\int{{{SC}_{s}^{g}\left( {t,a} \right)}{\delta\left( {{t - {{\hat{t}}_{s}\left( {t,a} \right)}},{a^{\prime} - {{\hat{a}}_{s}\left( {t,a} \right)}}} \right)}\frac{a^{\prime 2}{\mathbb{d}t}{\mathbb{d}a}}{a^{2}}}}}} & (22)\end{matrix}$benefits both from the smoothing performed by the mother wavelet, andfrom the reassignment, which refocuses the scalogram on the squeezedsignal description given by the Wigner-Ville distribution.

From a computational point of view, the local centroids can be computedefficiently by means of two additional wavelet transforms, using twoparticular mother wavelets:

$\begin{matrix}{{{{\hat{t}}_{s}\left( {t,a} \right)} = {t + {{Re}\left\{ \frac{{aW}_{s}^{tg}\left( {t,a} \right)}{W_{s}^{g}\left( {t,a} \right)} \right\}}}}{{{\hat{\omega}}_{s}\left( {t,a} \right)} = {\frac{\omega_{o}}{{\hat{a}}_{s}\left( {t,a} \right)} = {{- {Im}}\left\{ \frac{W_{s}^{d\; g}\left( {t,a} \right)}{{aW}_{s}^{g}\left( {t,a} \right)} \right\}}}}} & (23)\end{matrix}$Several mother wavelet functions can be used. One of them is the Morletwavelet also described previously and defined as follows:

$\begin{matrix}{{g(t)} = {\frac{1}{\sqrt{\sigma}}{\mathbb{e}}^{\frac{- t^{2}}{2\;\sigma^{2}}}{\mathbb{e}}^{{\mathbb{i}}\;\omega_{o}t}}} & (24)\end{matrix}$In that case, where dg(t) represents the derivative of the motherwavelet with respect to time,

$\begin{matrix}{{{{dg}(t)} = {\frac{- {{tg}(t)}}{\sigma^{2}} + {{\mathbb{i}}\;\omega_{o}{g(t)}\mspace{14mu}{and}}}}{{W_{s}^{dg}\left( {t,a} \right)} = {- \left\lbrack {\frac{W_{s}^{tg}\left( {t,a} \right)}{\sigma^{2}} + {{\mathbb{i}}\;\omega_{o}{W_{s}^{g}\left( {t,a} \right)}}} \right\rbrack}}} & (25)\end{matrix}$Only W_(s) ^(g)(t,a) and aW_(s) ^(tg)(t,a) need to be computed, becauseusing the expression in equation (21) yields:

$\begin{matrix}{{{\hat{\omega}}_{s}\left( {t,a} \right)} = {\frac{\omega_{o}}{{\hat{a}}_{s}\left( {t,a} \right)} = {\frac{\omega_{o}}{a} + {\frac{1}{a^{2}\sigma^{2}}{Im}\left\{ \frac{{aW}_{s}^{tg}\left( {t,a} \right)}{W_{s}^{g}\left( {t,a} \right)} \right\}}}}} & (26)\end{matrix}$Computer algorithms can therefore be generated from the discrete-timeversions of the following expressions:

$\begin{matrix}{{{W_{s}^{g}\left( {t,a} \right)} = {\sqrt{\frac{\omega }{\omega_{o}\sigma}}{\int_{- \infty}^{+ \infty}{{s\left( {t + \tau} \right)}{\mathbb{e}}^{\frac{{- \omega^{2}}\tau^{2}}{2\;\omega_{o}^{2}\sigma^{2}}}{\mathbb{e}}^{{- {\mathbb{i}}}\;\omega\;\tau}\ {\mathbb{d}\tau}}}}}{{{aW}_{s}^{tg}\left( {t,a} \right)} = {\sqrt{\frac{\omega }{\omega_{o}\sigma}}{\int_{- \infty}^{+ \infty}{\tau\;{s\left( {t + \tau} \right)}{\mathbb{e}}^{\frac{{- \omega^{2}}\tau^{2}}{2\;\omega_{o}^{2}\sigma^{2}}}{\mathbb{e}}^{{- {\mathbb{i}}}\;\omega\;\tau}\ {\mathbb{d}\tau}}}}}} & (27)\end{matrix}$with a=ω₀/ω. Since the Gaussian analyzing window used in theseexpressions depends on the frequency (or scale) parameter, FFTalgorithms cannot be used, resulting in much slower algorithms. Itshould be emphasized that ω_(o) and σ only appear through their productω_(o)σ, which is the only degree of freedom of this representation.Increasing this parameter improves the frequency resolution and reducesthe time resolution. The techniques described herein can be applied inthe same manner on the reassigned scalogram.

Dispersion Extraction From Wavelet Domain Problem Statement

It is useful to estimate the group and phase velocity of propagatingcomponents of acoustic data collected by an array of sensors. Suchacoustic data could be acquired in a variety of ways, including but notlimited to borehole acoustic tools, both wireline and LWD, and seismictools. In the illustrated example, the only assumption made about thedata is that it includes the superposition of one or more propagatingcomponents along with noise. Each of these components could exhibit bothattenuation and dispersion, and moreover could overlap with othercomponents in time or frequency or both.

The approach described here is based on the use of the complexcontinuous wavelet transform (CWT), such as described in the previoussection. The wavelet transform maps a one dimensional signal into a twodimensional (time-frequency) plane of coefficients. This increase indimensionality allows a better analysis and interpretation of the signaland a better separation of the various components. Thus, with realphysical signals, components are often resolved on the time-frequencyplane and can be separately reconstructed even when they do not separatein time or frequency. For the same reason, it is often possible toeffectively analyze the received signals using a single mode approach byoperating in the time-frequency or time-scale domain. In the embodimentdescribed below there is one propagating component exhibiting dispersiveand attenuative behavior. An example of an attenuative component wouldbe a leaky borehole mode, i.e., one that radiates into the formation.

For two sensors on an array, the Fourier transform of the signalreceived at the second sensor in terms of that at the first sensor canbe expressed as follows:

(f)=

(f)e ^(−2iπδk(f)) e ^(−δA(f)),   (28)where k(f) and A(f) are respectively the wavenumber and the attenuationas functions of frequency and δ is the inter-sensor spacing. Note thatin this case it is assumed that the signal s₁ arrives first (this is nota limitation however; if signal s₂ arrives first, equation (28) stillholds but with negative values for k(f) and A(f)). It is possible toexpress equation (28) linking the signals 1 and 2 in the continuouswavelet domain. The continuous wavelet transform can be expressed inFourier domain as follows:

$\begin{matrix}{{{S\left( {a,b} \right)} = {\int_{- \infty}^{+ \infty}{{\overset{\Cap}{g^{*}}({af})}{\mathbb{e}}^{2\;{\mathbb{i}}\;\pi\;{bf}}{\overset{\Cap}{s}(f)}\ {\mathbb{d}f}}}},} & (29)\end{matrix}$with

the complex conjugate of the Fourier transform of the analyzing motherwavelet and

the Fourier transform of the signal to analyze. Note that thenormalization coefficient has not been included for simplicity. Applyingthe wavelet transform to the equation (28) gives

$\begin{matrix}{{S_{2}\left( {a,b} \right)} = {\int_{- \infty}^{+ \infty}{{\overset{\Cap}{g^{*}}({af})}{\mathbb{e}}^{2\;{\mathbb{i}}\;\pi\;{bf}}{{\overset{\Cap}{s}}_{1}(f)}{\mathbb{e}}^{{- \mspace{2mu} 2}\;{\mathbb{i}}\;\pi\;\delta\;{k{(f)}}}\ {\mathbb{e}}^{{- \delta}\;{A{(f)}}}{\mathbb{d}f}}}} & (30)\end{matrix}$which can be simplified as

$\begin{matrix}{{S_{2}\left( {a,b} \right)} = {\int_{- \infty}^{+ \infty}{{\overset{\Cap}{g^{*}}({af})}{{\overset{\Cap}{s}}_{1}(f)}{\mathbb{e}}^{2\;{\mathbb{i}}\;\pi\;{f{\lbrack{b - \frac{\delta\;{k{(f)}}}{f}}\rbrack}}}{\mathbb{e}}^{{- \delta}\;{A{(f)}}}{\mathbb{d}f}}}} & (31)\end{matrix}$

In a more general way, it is possible to rewrite the above expression asa function of the sensor number, which yields:

$\begin{matrix}{{{S_{l}\left( {a,b} \right)} = {\int_{- \infty}^{+ \infty}{{\overset{\Cap}{g^{*}}({af})}{{\overset{\Cap}{s}}_{j}(f)}{\mathbb{e}}^{2\;{\mathbb{i}}\;\pi\;{f{\lbrack{b - \frac{\delta_{l,j}{k{(f)}}}{f}}\rbrack}}}{\mathbb{e}}^{{- \delta_{l,j}}{A{(f)}}}{\mathbb{d}f}}}},{l > j}} & (32)\end{matrix}$where δ_(l,j) now denotes the spacing between the l^(th) and j^(th)sensors. This equation gives the relationship between the signal l and jin the wavelet domain. The objective is now to extract the wavenumberand the attenuation in order to be able to derive the group and phasevelocity. There are two approaches to extract these parameters fromequation (32). The first approach explicitly inverts the equation usinga non-linear optimization approach, while the second approach, followedhere, involves taking a local linear expansion of the wavenumber andattenuation and solving the problem with the resulting simpler model. Inparticular, a local Taylor expansion of the parameters k (f) and A(f) isdeveloped. Numerical studies indicate that a first or second order localexpansion of the wavenumber dispersion is sufficient to capture thelocal variation; physical considerations support this by imposingsmoothness on the dispersion curves. FIGS. 2 a and 2 b illustrate theresults of approximating a typical dispersion curve (quadrupole) withfirst and second order Taylor expansions. The local fit obtained therebyis adequate for capturing the local behavior of the dispersion curve.Moreover, since the approximation used in equation (33) below needs tobe valid only over the spectral support of the scaled wavelet {tildeover (g)}(af) and since this interval is smaller at the scalescorresponding to low frequencies where the sharpest slowness variationtypically occurs, it usually suffices to use the linear (first order)Taylor expansion for the wavenumber dispersion.

Linear Approach

It is possible to express the local linear Taylor expansion of k(f) andA(f) as

$\begin{matrix}{{{k(f)} = {{k\left( f_{a} \right)} + {\left( {f - f_{a}} \right){k^{\prime}\left( f_{a} \right)}} + {O\left( {{f - f_{a}}} \right)}}}{and}} & (33) \\{{{A(f)} = {{A\left( f_{a} \right)} + {O\left( {{f - f_{a}}} \right)}}}{with}} & (34) \\{f_{a} = \frac{\omega_{0}}{2\;\pi\; a}} & (35)\end{matrix}$where a is the scale considered (i.e. the frequency) and ω₀ the centralfrequency of the mother wavelet.

Substituting the expressions (33) and (34) in (32) for the case of thelocal linear (order one) expansion yields:S _(t)(a,b)=e ^(−δ) ^(l,j) ^(A(f) ^(a) ⁾ e ^(−i2πδ) ^(i,j) ^([k(f) ^(a)^()−f) ^(a) ^(k′(f) ^(a) ^()]) S _(j)(a,b−δ _(l,j) k′(f _(a))).   (36)

Observing that

$\begin{matrix}\begin{matrix}{\left\lbrack {{k\left( f_{a} \right)} - {f_{a}{k^{\prime}\left( f_{a} \right)}}} \right\rbrack = {f_{a}\left\lbrack {\frac{k\left( f_{a} \right)}{f_{a}} - {k^{\prime}\left( f_{a} \right)}} \right\rbrack}} \\{= {{f_{a}\left\lbrack {{s_{\phi}\left( f_{a} \right)} - {s_{g}\left( f_{a} \right)}} \right\rbrack}\bullet\frac{\;\varphi_{a}}{2\;\pi}}}\end{matrix} & (37)\end{matrix}$appears as a phase factor, φ_(a), that can be written in terms of thephase and group slowness; taking j=l₀, the reference sensor, andrepresenting the distance to it of the l^(th) sensor by δ_(l), andfinally writing the time shift parameter b simply as time, t, in otherwords, considering the wavelet coefficients as waveforms in time, it ispossible to rewrite equation (36) asS _(l)(a,t)=e ^(−δ) ^(l) ^(A(f) ^(a) ⁾ e ^(−iδ) ^(l) ^(φ) ^(a) S _(l) ₀(a,t−δ _(l) s _(g))   (38)

Equation (38) shows the relationship between the wavelet transformcoefficients at each scale a of a pair of signals (i.e. l and l₀) usingthe linear Taylor development of the wavenumber dispersion. In that caseit will be recognized that the CWT coefficients at the l^(th) sensor aretime shifted with respect to those at the l₀ ^(th) sensor by an amountgiven by the group slowness and inter-sensor distance, and multiplied bya factor whose magnitude is dependent on the attenuation and whose phasedepends on the difference of the phase and group slowness. Therefore,given the coefficients at a particular scale, a, corresponding to thefrequency, f_(a), there are three quantities to estimate, namely, theattenuation factor, the phase factor and the time shift given by thegroup slowness.

Method-1: Extracting the Group Slowness from the Modulus Maxima of theWavelet Transform

Group slowness represents the velocity with which the wave's envelope(shape of the amplitude) and energy propagates through space. In thecase of a dispersive or attenuating medium this group velocity becomesdispersive, i.e., a function of the frequency. As already described, thetransformation conserves the energy of the signal and therefore it ispossible to obtain an estimate of the group velocity as a function offrequency directly in the wavelet domain. It will now be demonstratedhow it is possible to extract the group velocity of a signal using thelocal modulus maxima of the wavelet transform.

In general, a wave propagating in a plane can be specified by thewavelength λ and the angle Θ between the x₁ axis and the propagationdirection. Hence, it is possible to obtain the apparent wavelength inthe x₁ and x₂ direction respectively, by

$\lambda_{1} = {{\frac{\lambda}{\cos\;\Theta}\mspace{14mu}{and}\mspace{14mu}\lambda_{2}} = \frac{\lambda}{\sin\;\Theta}}$(See FIG. 3). The superposition of two waves traveling in a plane withthe same unit amplitude and with frequencies ω₁ and ω₂, respectively, isgiven byu(x ₁ , x _(2,) t)=e ^(i(ω) ¹ ^(t−k) ¹¹ ^(x) ¹ ^(−k) ¹² ^(x) ² ⁾ +e^(i(ω) ² ^(t−k) ²¹ ^(x) ¹ ^(−k) ²² ^(x) ² ⁾   (39)where k_(ij) are wavenumbers corresponding to the frequency ω_(j) and tothe coordinate x_(i) with the wavenumbers defined as

$k = {\frac{2\;\pi}{\lambda}.}$Applying this definition to waves traveling in the same direction leadsto the wavenumbers k_(j1)=k_(j) cos Θ and k_(j2)=k_(j) sin Θ. Thereforethe wave equation can be expressed asu(x ₁ ,x _(2,) t)=e ^(i(ω) ¹ ^(t−k) ¹ ^(cos Θx) ¹ ^(−k) ¹ ^(sin Θx) ² ⁾+e ^(i(ω) ² ^(t−k) ² ^(cos Θx) ¹ ^(−k) ² ^(sin Θx) ² ⁾   (40)Defining

${k_{c} = \frac{k_{1} + k_{2}}{2}},{\omega_{c} = \frac{\omega_{1} + \omega_{2}}{2}},{{\Delta\; k} = \frac{k_{1} - k_{2}}{2}},{{{and}\mspace{14mu}\Delta\;\omega_{c}} = \frac{\omega_{1} - \omega_{2}}{2}},$the wave packet becomesu(x ₁ ,x ₂ ,t)=2 cos(Δω_(t) −Δk cos Θx ₁ −Δk sin Θx ₂)e ^(i(ω) ^(c)^(t−k) ^(c) ^(cos Θx) ¹ ^(−k) ^(c) ^(sin Θx) ² ⁾   (41)Fourier and wavelet transform of u(x₁,x₂,t) yields:U(x ₁ ,x ₂ ,a,b)=√{square root over (a)}[e ^(i(ω) ¹ ^(b−k) ¹ ^(cos Θx) ¹^(−k) ¹ ^(sin Θx) ² ^() g(ω) ¹ ⁾ +e ^(i(ω) ² ^(b−k) ² ^(cos Θx) ¹ ^(−k)² ^(sin Θx) ² ^() g) *^((aω) ¹ ⁾]  (42)After manipulation, the result is:

$\begin{matrix}{{{U\left( {x_{1},x_{2},a,b} \right)}} = {\sqrt{a}\begin{bmatrix}{{\overset{\Cap}{g}\left( {a\;\omega_{1}} \right)}^{2} + {\overset{\Cap}{g}\left( {a\;\omega_{2}} \right)}^{2} + {2\;{\overset{\Cap}{g}\left( {a\;\omega_{1}} \right)}}} \\{{\overset{\Cap}{g}\left( {a\;\omega_{2}} \right)}{\cos\left( {{2\;\Delta\;\omega_{b}} - {2\;\Delta\; k\;\cos\;\Theta\; x_{1}} - {2\;\Delta\; k\;\sin\;\Theta\; x_{2}}} \right)}}\end{bmatrix}}^{\frac{1}{2}}} & (43)\end{matrix}$If k is small, then k₁≈k₂, ω₁≈ω₂ and ĝ(aω₁)≈ĝ(aω₂)≈ĝ(aω_(c)), then (43)becomes

$\begin{matrix}{{{U\left( {x_{1},x_{2},a,b} \right)}} \approx {\sqrt{2\; a}{{{{\overset{\Cap}{g}\left( {a\;\omega_{c}} \right)}}\left\lbrack {1 + {2\;{\cos\left( {{2\;\Delta\;\omega_{b}} - {2\;\Delta\; k\;\cos\;\Theta\; x_{1}} - {2\;\Delta\; k\;\sin\;\Theta\; x_{2}}} \right)}}} \right\rbrack}^{\frac{1}{2}}.}}} & (44)\end{matrix}$If

$\begin{matrix}{{{{2\;\Delta\;\omega_{b}} - {2\;\Delta\; k\;\cos\;\Theta\; x_{1}} - {2\;\Delta\; k\;\sin\;\Theta\; x_{2}}} = 0}{and}} & (45) \\{a = \frac{\omega_{0}}{\omega_{c}}} & (46)\end{matrix}$the magnitude of the wavelet transform takes its maximum value. This isthe case for

$\begin{matrix}{b = {{\frac{\Delta\; k}{\Delta\;\omega}\left( {{\cos\;\Theta\; x_{1}} + {\sin\;\Theta\; x_{2}}} \right)} = {\frac{1}{c}{\left( {{\cos\;\Theta\; x_{1}} + {\sin\;\Theta\; x_{2}}} \right).}}}} & (47)\end{matrix}$

This means that the location of the maximum peak of the magnitude of thewavelet transform at scale a provides the arrival time of thepropagating wave with a group velocity c at the corresponding frequency.Therefore the maxima line of the wavelet transform will give the arrivaltime of the component of interest as a function of the frequency, whichwill provide in a straightforward manner the group velocity averagedover the source to sensor distance. In order to extract the groupslowness across the array of sensors, it is possible to fit a line inthe least squares sense to the arrival times at the sensors for eachscale and use the slope of the fitted line to obtain the group slownessestimate at the corresponding center frequency, f_(a). The peaklocations can be corrected via exponential quadratic interpolationacross time. The exponential quadratic interpolation may be chosen sincethe envelope of the reproducing kernel is a quadratic exponential in‘b’, (Grossmann, A. and Morlet, J., 1984, Decomposition of Hardyfunctions into square integrable wavelets of constant shape, SIAM—J.Math. Anal., 15, 723-736).

FIGS. 4, 5 and 6 illustrate a methodology to extract the group slownessfrom the modulus of the wavelet transform of the recorded waveforms.

From Group Slowness to Phase Slowness and Attenuation

In order to obtain phase slowness from group slowness, a shift is firstapplied to the wavelet coefficients at each frequency, given by(δ_(l)s_(g)(f_(a))) as per equation (38), using the estimates of thegroup slowness obtained as described above. After the shift is applied,one approach for obtaining the phase factor, φ_(a), is to take theFourier transform along the sensor axis on the shifted waveletcoefficients at the location of the intercept of the fitted line on thereference sensor or in a window around that location. The maximum of theFFT peak will correspond to the phase factor φ_(a). Using therelationship presented previously i.e. φ_(a)=2πf_(a)[s₁₀₀−s_(g)], it ispossible to obtain the phase slowness.

The FFT operation can be viewed as a fast search to find the desiredphase correction and is applicable for arrays with uniformly spacedsensors on a line, also called a uniform linear array or ULA. Theattenuation too can be found by a similar search operation over possiblevalues of the attenuation, but this can become computationallyintensive. An alternative approach which simultaneously estimates thegroup and phase slowness is described below.

Assuming as above that the correct time shift is applied to the waveletcoefficients, S_(l)(a,t), a rank one subspace model is obtained for theshifted coefficients:

$\begin{matrix}{\begin{bmatrix}{S_{1}\left( {a,{t^{\prime} + {\delta_{1}s_{g}}}} \right)} \\{S_{2}\left( {a,{t^{\prime} + {\delta_{2}s_{g}}}} \right)} \\\vdots \\{S_{L}\left( {a,{t^{\prime} + {\delta_{M}s_{g}}}} \right)}\end{bmatrix} = {\begin{bmatrix}{\mathbb{e}}^{- {\delta_{1}{({{A{(f_{a})}} + {{\mathbb{i}}\;\varphi_{a}}})}}} \\{\mathbb{e}}^{- {\delta_{2}{({{A{(f_{a})}} + {{\mathbb{i}}\;\varphi_{a}}})}}} \\\vdots \\{\mathbb{e}}^{- {\delta_{L}{({{A{(f_{a})}} + {{\mathbb{i}}\;\varphi_{a}}})}}}\end{bmatrix}{S_{l_{0}}\left( {a,t^{\prime}} \right)}}} & (48)\end{matrix}$where l₀ refers to the reference sensor and t′ refers to a set of timeindices in a window encompassing the mode of interest.

Note that this formulation corresponds to the complex exponentialestimation problem. While, in general, this problem requires anon-linear inversion (or 2-D search), faster methods have been developedfor the special but common case of the uniform linear array (ULA). Forexample, the matrix pencil method, such as in Ekstrom, M. P., 1995,Dispersion estimation from borehole acoustic arrays using a modifiedmatrix pencil algorithm, volume 2 of 29th Asilomar Conference onSignals, Systems and Computers, pages 449-453, could be used to estimatethe exponential parameter.

Alternatively, an even faster approach may be utilized which is easilyapplied to the ULA (and also to the case when the array can besubdivided into two sub-arrays with constant inter-sensor spacingbetween the corresponding elements of each), and which yields a closedform solution. However, it should be appreciated that it is notnecessary to have a ULA for this method.

For purposes of this explanation, δ is defined as the commoninter-sensor spacing between adjacent sensors, noting thatδ_(l)=δ(l−l₀), dropping the subscript on φ=2πf_(a)(s_(φ)−s_(g)), andfurther defining α=A(f_(a)), and Y_(a) as the LHS of equation (48), andrewrite that as:

$\begin{matrix}{Y_{a} = {{\begin{bmatrix}{\mathbb{e}}^{{- {\delta{({\alpha + {{\mathbb{i}}\; 2\;\pi\;\varphi}})}}}{({1 - l_{0}})}} \\{\mathbb{e}}^{{- {\delta{({\alpha + {{\mathbb{i}}\; 2\;\pi\;\varphi}})}}}{({2 - l_{0}})}} \\\vdots \\{\mathbb{e}}^{{- {\delta{({\alpha + {{\mathbb{i}}\; 2\;\pi\;\varphi}})}}}{({L - l_{0}})}}\end{bmatrix}{S_{l_{0}}\left( {a,t^{\prime}} \right)}} + N}} & (49)\end{matrix}$(considering the presence of noise, N).

Defining

$\begin{matrix}{{Y_{a,1} = \begin{bmatrix}{Y_{a}(1)} \\{Y_{a}(2)} \\\vdots \\{Y_{a}\left( {L - 1} \right)}\end{bmatrix}},{Y_{a,2} = \left\lbrack \begin{matrix}{Y_{a}(2)} \\{Y_{a}(3)} \\\vdots \\{Y_{a}(L)}\end{matrix} \right)}} & (50)\end{matrix}$compute the quantities

$\begin{matrix}{R_{ij} = {\sum\limits_{o}{Y_{a,i}^{*}\bullet\; Y_{a,j}}}} & (51)\end{matrix}$for i,j=1,2 where (·)* denotes the complex conjugate, □ implies theelement-by-element product of the matrices, Y_(a,i) and Y_(a,j) andΣ_(o) indicates a sum taken over the elements of the product matrix soobtained. Note that R₁₂ and R₂₁ are complex conjugate and so only one ofthose needs to be computed in practice.

Provided that the noise is white Gaussian (which is reasonable if thereis no coherent interference at the chosen moveout), the matrix built upfrom the above quantities satisfies

$\begin{matrix}{\begin{bmatrix}R_{11} & R_{12} \\R_{21} & R_{22}\end{bmatrix} = {{\sigma_{1}V^{H}} + {\sigma_{2}I_{2 \times 2}}}} & (52)\end{matrix}$where V=[1 e^(−δ(α+iφ))], I_(2×2) is the identity matrix of size 2 andσ₁ and σ₂ are two positive constants. Therefore, V is the righteigenvector corresponding to the larger (dominant) eigenvalue. Also, itscomponents can be obtained in closed form for this case. Results for theestimates of α and φ can be expressed as:

$\begin{matrix}{{\overset{\Cap}{\alpha} = {{- {\log\left( \frac{\sqrt{\left( {R_{11} - R_{22}} \right)^{2} + {4{R_{12}}^{2}}} - R_{11} + R_{22}}{2{R_{12}}} \right)}}/\delta}}{\overset{\Cap}{\varphi} = {{- {\angle\left( R_{12} \right)}}/\delta}}} & (53)\end{matrix}$where ∠(R₁₂) denotes the complex phase of R₁₂.

Note that while the above process is designed to be applied to a ULA oran array that can be decomposed into two sub-arrays with constantseparation, it is possible to use a similar process for a general arrayto estimate the attenuation and phase factors. In particular, theprocess is applied to the sensors in pairwise fashion, using theappropriate inter-sensor spacing δ for each pair and taking the average.

As already explained above, these estimates can now be used to generateestimates for the phase slowness and attenuation at the frequencycorresponding to the scale a being processed. Repeating this for theother scales (frequencies) of interest yields the desired dispersioncurve estimates.

Multiple Peaks Processing

As outlined above, in order to find the group slowness at eachfrequency, the shift in the peak at the particular frequency isdetermined, i.e., at a particular value of the scale a. In the generalset-up there can be more than one ‘mode’ that can be present. Thus, ingeneral:

$\begin{matrix}{{S_{l}\left( {a,\; b} \right)} = {\sum\limits_{i = 1}^{P}{S_{l}^{i}\left( {a,b} \right)}}} & (54)\end{matrix}$Assuming that the time frequency separation of the modes is sufficientto allow the corresponding energy concentrations to occupy disjointregions (or support sets) in the time-frequency plane across the arrays,the following approximation holds true:

$\begin{matrix}{{{S_{l}\left( {a,b} \right)}} \sim {\sum\limits_{i = 1}^{P}{{S_{l_{0}}^{i}\left( {a,{b - {{s_{g}^{i}\left( f_{a} \right)}\delta_{l}}}} \right)}}}} & (55)\end{matrix}$In particular, equation (55) indicates that it is possible to obtain theenvelope peaks for each of the modes when the modes are sufficientlyseparated in the time frequency place, for accurate processing of groupslowness for each mode, in presence of low noise.

Data Association

In order to associate peaks from one sensor to the next, it is assumedthat the group slowness is not varying across the sensors. Since themodes are well separated in time and frequency (i.e., the peaks are wellseparated), the problem is formulated as a 2-D minimal matching problemfrom a current set of peaks to the next set of peaks. This can beaccomplished as follows:

Fix a frequency, e.g., f_(a). Let P be the number of peaks that arebeing tracked at frequency f_(a). Consider a sensor l and denote thepeak locations at sensor l by Y_(i),i=1, 2, . . . , P. Denote the peaklocations at sensor l+1 by Y_(j), j=1, 2, . . . , P. Then, to solve thefollowing minimization problem,min_(e) _(ij) Σ_(i, j)c_(ij)e_(ij)sub to: e_(ij) ∈0, 1; Σ_(i)e_(i,j)=1; Σ_(j)e_(ij)=1   (56)where e_(ij)=0 denotes no association between Y_(i) and Y_(j) ande_(ij)=1 denotes that Y_(i) is associated with Y_(j). The cost c_(ij) isproportional to the dissimilarity between Y_(i) and Y_(j) and can becaptured via many metrics. It is also possible to incorporate ‘prior’information in order to develop the cost structure. The above integerprogramming problem can be relaxed to a linear programming one via thefollowing approach. Let A denote the vertex edge incidence matrix forthe complete bipartite graph between vertex sets Y_(i) and Y_(j), i.e.A_(ve)=1 if vertex v ∈ edge(e), otherwise A_(ve)=0. With the notation1=[1, 1, . . . , 1]^(t) and e=[e₁₁e₁₂, . . . , e_(ik), . . . ]^(t), itis possible to write the above integer program asmin c^(t)esub to:Ae=1; e≧0   (57)where c is a N²×1 vector of costs on edge e_(ij), ∀i,j. The theoremdescribed below is useful for this purpose. (Steele, J. M., Probabilitytheory and combinatorial optimization, SIAM).Theorem

For the vertex edge incidence matrix A of a bipartite graph, thevertices of the polytope defined byAe≦1; e≧0consists of 0-1 vectors. Moreover these vectors are incidence vectors ofmatchings.□

From the above theorem, it follows that the integer programming problemcan be relaxed to a linear programming one by relaxing the hardconstraints of equalities to inequalities. Then using the “SIMPLEX”method (O(P³) complexity) of linear programming it is possible to findthe minimal cost matching.

In the method described herein it is assumed that the group slownessdoes not vary significantly across the sensor array, that it is boundedfrom above, and there is no negative group velocity. Based on theseassumptions, the costs are chosen as

$\begin{matrix}{c_{ij} = \left\{ \begin{matrix}{Y_{j} - Y_{i}} & {{{{if}\mspace{14mu} Y_{j}} - Y_{i}} > 0} \\C_{0} & {{{{if}\mspace{14mu} Y_{j}} - Y_{i}} < {{0\mspace{14mu}{or}\mspace{14mu}{if}\mspace{14mu} Y_{j}} - Y_{i}} > T_{0}}\end{matrix} \right.} & (58)\end{matrix}$where C₀ is a large number and T₀ is the limit on the time shift and canbe derived from the sensor spacing and the sampling time and upper boundon the group slowness. Other choices of cost function are also possibleand several lead to good results in practice.

Propagating the Association Across Sensors

Associating the peaks across the sensors, 1, . . . , L is an associationproblem over L sets Y₁, . . . , Y_(L) of peak locations. This is aL-dimensional data association problem which can be formulated as aninteger program in a similar manner. However the complexity of such analgorithm is o((P!)^(L)).

Suboptimal Approach

Since the complexity of association across the sensors is exponential,it may be preferable to employ a suboptimal approach for dataassociation across the sensors. One example of a suboptimal algorithm isas follows.

Algorithm

1. Start from the last sensor as the modes have more chance to beseparated in time at the last sensor. Let the peak locations be given byY_(L). Assume these to be the correct peaks.

2. Using 2-D minimal assignment, find the correct permutation, π_(L−1),for the peak locations, Y_(L−1), with respect to the peak locations,Y_(L). Call the new permuted vector Y_(L−1)(π_(L−1)).

3. Using 2-D minimal assignment, find π_(L−2) for peaks Y_(L−2) withrespect to Y_(L−1)(π_(L−1)).

4. Repeat steps 2 and 3 until the first sensor.

The above algorithm behaves like an inference propagation type algorithmwhere the inference from the previous association is propagated to thenext. This approach is suboptimal because it is subject to propagationof errors. However, error propagation can be mitigated by using goodcost structures, i.e., incorporating some prior information or if thenumber of peaks are small then using K dimensional matching and thenpropagating the associations.

The described data association process can also be applied to trackingpeaks across frequency (or scale) on the time frequency plane. This isuseful for labeling propagating modes when more than one are present,and using that to appropriately extract the dispersion curve for each.FIG. 7 illustrates an example of such a mode tracking on the timefrequency plane using the data association described above.

Bias Correction

Estimates of the dispersion can be biased in some cases. The analysisfor the method described below shows that a bias may exist particularlyfor the group slowness when the spectral mean frequency of the waveletcoefficients at any scale are different from the nominal centerfrequency and when the group slowness changes over the bandwidth of thesaid coefficients. A relatively simple technique to correct theresulting bias is to refer the dispersion estimates to the spectral meanfrequency rather than to the nominal center frequency. The spectral meanis calculated as follows

$\begin{matrix}{f_{c} = \frac{\int{\sum\limits_{l}{{{{\overset{\sim}{S}}_{l}\left( {a,f} \right)}}^{2}f{\mathbb{d}f}}}}{\int{\sum\limits_{l}{{{{\overset{\sim}{S}}_{l}\left( {a,f} \right)}}^{2}{\mathbb{d}f}}}}} & (59)\end{matrix}$where S(a,f) is the Fourier transform of the wavelet coefficientscalculated in a neighborhood of the mode of interest. This could also beobtained as the instantaneous frequency computed as a derivative of thephase of the complex coefficients at the modulus maxima. Note that whenthe coefficients from the reassigned scalogram are used, this frequencycorrection is already applied as a consequence of the reassignment asshown in equation (21) and no further bias correction is generallyneeded.

FIG. 8 illustrates an example of application of the proposed method toreal oilfield data. In both cases, the proposed method is compared tothe existing matrix pencil method. Note that the new method yields morestable and better estimates.

Method-2: The Exponential Projected Radon Transform (EPRT)

As indicated by equation (38), the three quantities to estimate are theattenuation factor, the phase factor and the time shift given by thegroup slowness. The last problem is similar to non-dispersive processingwhich has been well studied in the literature. A classic method is touse the Radon transform (Deans, S. R., 1983, The Radon Transform andSome of its Applications, A Wiley Interscience publication) or operatorsbased thereon. The Radon transform is a way of representing an array ofwaveform traces acquired at a number of known spatial locations in termsof components propagating at different moveouts or slownesses. This canbe used to separate non-dispersive components propagating at differentslownesses and in addition can also be used to estimate thesequantities. Thus, if x_(l)(t) represents a set of waveform tracesobtained with an array of sensors and consisting of non-dispersivecomponents satisfying the relationy _(l) ^(i)(t)=y _(l) _(o) ^(i)(t−s _(i)δ_(l))   (60)where y^(i) is the i^(th) non-dispersive component of x propagating atthe slowness s_(i) and the Radon transform takes the form of a slantstack:

$\begin{matrix}{{R\left( {p,t} \right)} = {\sum\limits_{l = 1}^{L}{x_{l}\left( {t + {p\;\delta_{l}}} \right)}}} & (61)\end{matrix}$where p is the test moveout or slowness.

Associated with this, it is possible to compute a semblance, which iscalculated by taking the ratio of the Radon transform output to thetotal energy. In practice, in order to account for the presence of morethan one component in the input data and to improve the noiseperformance, this ratio is taken after integrating both the slant stackenergy (which captures the energy of the coherent component at thecorresponding moveout) and the input energy over a set of time windowschosen to capture the individual components and indexed by the startingtime

$\begin{matrix}{{\rho\left( {t,p} \right)} = \frac{\int_{t}^{t + T_{w}}{{{\sum\limits_{l = 1}^{L}{x_{l}\left( {t + {p\;\delta_{l}}} \right)}}}^{2}\ {\mathbb{d}t}}}{L{\int_{t}^{t + T_{w}}{\sum\limits_{l = 1}^{L}{{{x_{l}\left( {t + {p\;\delta_{l}}} \right)}}^{2}\ {\mathbb{d}t}}}}}} & (62)\end{matrix}$The semblance can therefore be interpreted as the ratio of the coherentenergy divided by the total energy in a time window T_(w) at time t andmeasures the degree of similarity of the received waveforms in thewindows corresponding to the test moveout across the array and hence thename. When there is a perfect match of the test moveout p with theslowness s_(i) of a component, the result is a value close to 1 which isthe maximum possible. In other words, the peaks of ρ in the (t, p) planeyield the slowness and time locations of compact non-dispersivepropagating components. This forms the basis of theslowness-time-coherence processing (Kimball, C. V. and Marzetta, T.,Semblance Processing of borehole acoustic array data, Geophysics,49(3):264-281, March 1984) applied to non-dispersive waves in theoilfield business.

The Exponential Projected Radon Transform (EPRT)

Although estimating phase and attenuation factors in addition to thesimple time shift of the non-dispersive processing above is morecomplex, it is possible to address the issue of estimating theseadditional parameters by noting that the Radon transform and semblancedescribed above perform the operation of projecting onto a subspaceafter performing a test moveout correction by time-shifting thewaveforms. In that case, the subspace is a special one comprising avector of ‘1’s. When the test moveout matches the true slowness, thenthe projection on this subspace leads to a maximum of the chosencriterion, either energy or semblance.

Equation (38) presents an analogous situation. Assuming that it ispossible to apply the correct timeshift corresponding to the groupslowness, s_(g), to the wavelet coefficients, S_(l)(a, t), then a rankone subspace model for the shifted coefficients can be obtained as shownpreviously in equation (48). Referring to that equation, by analogy itis desirable to project onto the subspace corresponding to the columnvector containing the exponential terms on the RHS. However, theattenuation and phase parameters comprising the latter are unknown andit is first necessary to estimate them. An approach to do so, given atime shift, has already been described under method 1 with the desiredestimates being given by equation (53).

Having thus obtained these exponential parameter estimates, it is nowpossible to perform the desired subspace projection as described above.Further, making the ULA assumption, it is possible to describe theestimated subspace (from equation 48) as:

$\begin{matrix}{U = \begin{bmatrix}{\mathbb{e}}^{{- {({\overset{\Cap}{\alpha}\; + {{\mathbb{i}}\;\overset{\Cap}{\varphi}}})}}{\delta{({1 - l_{0}})}}} \\{\mathbb{e}}^{{- {({\overset{\Cap}{\alpha}\; + {{\mathbb{i}}\;\overset{\Cap}{\varphi}}})}}{\delta{({2 - l_{0}})}}} \\\vdots \\{\mathbb{e}}^{{- {({\overset{\Cap}{\alpha}\; + {{\mathbb{i}}\;\overset{\Cap}{\varphi}}})}}{\delta{({L - l_{0}})}}}\end{bmatrix}} & (63)\end{matrix}$The projection operator on this subspace is give by:

$\begin{matrix}\begin{matrix}{P_{U} = {\frac{1}{\sqrt{U^{H}U}}U^{H}}} \\{= {\frac{1}{\sqrt{\sum\limits_{l}{\mathbb{e}}^{{- 2}\overset{\Cap}{\;\alpha}\;{\delta{({l - l_{0}})}}}}}U^{H}}} \\{= {\left( {{\mathbb{e}}^{{- \overset{\Cap}{\alpha}}\;{\delta{({{2\; l_{0}} + 1 - L})}}}\frac{\sinh\left( {\overset{\Cap}{\alpha}\;\delta} \right)}{\sinh\left( {L\;\overset{\Cap}{\alpha}\;\delta}\; \right)}} \right)^{1/2}U^{H}}}\end{matrix} & (64)\end{matrix}$where (□)^(H) denotes the Hermitian or complex conjugate transpose.

Applying this projection on the wavelet coefficients computed on thearray data at a particular scale a for a set of time and moveouts yieldsthe analogous form for the modified Radon transform

$\begin{matrix}{{R_{a}\left( {t,{p;\overset{\Cap}{\alpha}},\overset{\Cap}{\varphi}} \right)} = {{\mathbb{e}}^{- {\overset{\Cap}{\alpha}{({{2\; l_{0}} + 1 - L})}}}\frac{\sinh\left( {\overset{\Cap}{\alpha}\;\delta} \right)}{\sinh\left( {L\;\overset{\Cap}{\alpha}\;\delta} \right)}{\int_{t}^{t + T_{w}}{{{\sum\limits_{l = 1}^{L}{{\mathbb{e}}^{{- {({\overset{\Cap}{\alpha} - {{\mathbb{i}}\;\overset{\Cap}{\varphi}}})}}{\delta{({l - l_{0}})}}}{S_{l}\left( {a,{t + {{p\left( {l - l_{0}} \right)}\delta}}} \right)}}}}^{2}\ {\mathbb{d}t}}}}} & (65)\end{matrix}$where the energy is considered and therefore the projected quantitiesare squared and this energy is integrated in a window positionedaccording to the parameter t. The quantities

and

are now estimated for each t and p and therefore functions of thelatter. This operation is illustrated in FIG. 9.

As before, it is possible to compute a corresponding semblance quantityfor each scale

$\begin{matrix}{{\rho_{a}\left( {t,{p;\overset{\Cap}{\alpha}},\overset{\Cap}{\varphi}} \right)} = {{\mathbb{e}}^{- {\overset{\Cap}{\alpha}{({{2\; l_{0}} + 1 - L})}}}\frac{\sinh\left( {\overset{\Cap}{\alpha}\;\delta} \right)}{\sinh\left( {L\;\overset{\Cap}{\alpha}\;\delta} \right)}\frac{\int_{t}^{t + T_{w}}{{{\sum\limits_{l = 1}^{L}{{\mathbb{e}}^{{- {({\overset{\Cap}{\alpha} - {{\mathbb{i}}\;\overset{\Cap}{\varphi}}})}}{\delta{({l - l_{0}})}}}{S_{l}\left( {a,{t + {{p\left( {l - l_{0}} \right)}\delta}}} \right)}}}}^{2}\ {\mathbb{d}t}}}{\int_{t}^{t + T_{w}}{\sum\limits_{l = 1}^{L}{{{S_{l}\left( {a,{t + {{p\left( {l - l_{0}} \right)}\delta}}} \right)}}^{2}\ {\mathbb{d}t}}}}}} & (66)\end{matrix}$where when

=0, the expression is recovered for the non-dispersive case. It is thuspossible to obtain maps on the (t, p) plane analogous to the Radontransform and semblance where each point on the map is computed usingthe corresponding estimated quantities

and

for the projection. These maps are referred to as the ExponentialProjected Radon Transform (EPRT) and the EPRT semblance. As before, thepeaks give information about the time localization and group slowness ofthe propagating components at the scale being analyzed while thecorresponding estimated phase and attenuation factors can be used toextract the corresponding phase slowness and attenuation. Theseestimates can be further refined, if desired, by executing a localsearch in the vicinity of these estimates and the corresponding peakvalues of t and p to find the phase and attenuation factors thatmaximize the output energy and semblance respectively.

This extraction and possible refinements are discussed further in thenext section where these outputs are analyzed.

The Analysis of EPRT

For the sake of simplicity it is assumed that there is only one modepresent. Following the analysis in C. V. Kimbal and T. MarzettaSemblance Processing of borehole acoustic array data., Geophysics,49(3):264-281, March 1984, assume that the time window length T_(w)(note that this depends on the scale) and time position t at scale a ischosen in equation (66) so that it encloses the coefficient waveform ofinterest completely. Then using Parseval's relation it is possible towrite

$\begin{matrix}{{\rho_{a}\left( {{p;\overset{\Cap}{\alpha}},\overset{\Cap}{\varphi}} \right)} = {{\mathbb{e}}^{- {\overset{\Cap}{\alpha}{({{2\; l_{0}} + 1 - L})}}}\frac{\sinh\left( {\overset{\Cap}{\alpha}\;\delta} \right)}{\sinh\left( {L\;\overset{\Cap}{\alpha}\;\delta} \right)}\frac{\int{{{\sum\limits_{l = 1}^{L}{{S_{l}\left( {a,f} \right)}{\mathbb{e}}^{{- {({\overset{\Cap}{\alpha} - {{\mathbb{i}}{({\overset{\Cap}{\varphi} + {2\;\pi\;{pf}}})}}})}}{\delta{({l - l_{0}})}}}}}}^{2}{\mathbb{d}f}}}{\int{\sum\limits_{l = l}^{L}\;{{{S_{l}\left( {a,f} \right)}}^{2}{\mathbb{d}f}}}}}} & (67)\end{matrix}$where S_(l)(a, f) is the Fourier Transform of S_(l)(a, t) and ρ_(a)(p;{circumflex over (α)}, {circumflex over (φ)}) is the value of thesemblance for a value of t, the time window position, so as to encompassthe signal appropriately. As before, let k(f) be the actual (real)wavenumber as a function of the frequency for the dispersive wave andA(f) be the attenuation. Then it follows from the linearity of thewavelet transform thatS _(l)(a, f)=X(a, f)e ^(−(A(f)+i2πk(f))δ(l−l) ⁰ ⁾   (68)where X(a, f)=X(f)G*(af), G(f) is the Fourier transform of the analyzingmother wavelet and X(f) is the Fourier transform of the waveformcomponent contained in the applied time window and received at thereference sensor. Inserting this expression into the semblance computedby the EPRT (equation 67) yields

$\begin{matrix}{{\rho_{a}\left( {{p;\overset{\Cap}{\alpha}},\overset{\Cap}{\varphi}} \right)} = {\frac{1}{K}\frac{\int{{{\sum\limits_{l = 1}^{L}{{X\left( {a,f} \right)}{\mathbb{e}}^{{- {({{A{(f)}} + {{\mathbb{i}}\; 2\;\pi\;{k{(f)}}}})}}{\delta{({l - l_{0}})}}}{\mathbb{e}}^{{- {({\overset{\Cap}{\alpha} - {{\mathbb{i}}{({\overset{\Cap}{\varphi} + {2\;\pi\;{pf}}})}}})}}{\delta{({l - l_{0}})}}}}}}^{2}{\mathbb{d}f}}}{\int{\sum\limits_{l = 1}^{L}{{{{X\left( {a,f} \right)}{\mathbb{e}}^{{- {({{A{(f)}} + {{\mathbb{i}}\; 2\;\pi\;{k{(f)}}}})}}{\delta{({l - l_{0}})}}}}}^{2}{\mathbb{d}f}}}}}} & (69)\end{matrix}$where the normalizing quantity K (function of

) is defined by

$\begin{matrix}{\frac{1}{K} = {{\mathbb{e}}^{- {\overset{\Cap}{\alpha}{({{2\; l_{0}} + 1 - L})}}}\frac{\sinh\left( {\overset{\Cap}{\alpha}\;\delta} \right)}{\sinh\left( {L\;\overset{\Cap}{\alpha}\;\delta} \right)}}} & (70)\end{matrix}$Collecting attenuation and phase terms in equation (69) and simplifyingyields

$\begin{matrix}{{\rho_{a}\left( {{p;\overset{\Cap}{\alpha}},\overset{\Cap}{\varphi}} \right)} = {\frac{1}{K}\frac{\int{{{X\left( {a,f} \right)}}^{2}{{\sum\limits_{l = 1}^{L}{\mathbb{e}}^{{\{{{- {({\overset{\Cap}{\alpha} + {A{(f)}}})}} + {{\mathbb{i}}{({\overset{\Cap}{\varphi} + {2\;{\pi{({{pf} - {k{(f)}}})}}}})}}}\}}{\delta{({l - l_{0}})}}}}}^{2}{\mathbb{d}f}}}{\int{{{X\left( {a,f} \right)}}^{2}{\sum\limits_{l = 1}^{L}{{{\mathbb{e}}^{{- {A{(f)}}}{\delta{({l - l_{0}})}}}}^{2}{\mathbb{d}f}}}}}}} & (71)\end{matrix}$Denoting

$\begin{matrix}{{\Pi(f)} = {{\sum\limits_{l = 1}^{L}{{\mathbb{e}}^{{- {A{(f)}}}{\delta{({l - l_{0}})}}}}^{2}} = {{\mathbb{e}}^{{- \delta}\;{A{(f)}}{({L - 1 - {2\; l_{0}}})}}\frac{\sinh\left( {L\;\delta\;{A(f)}} \right)}{\sinh\left( {\delta\;{A(f)}} \right)}}}} & (72) \\{and} & \; \\{{\overset{\sim}{X}\left( {a,f} \right)} = {{X\left( {a,f} \right)}\sqrt{\Pi(f)}}} & (73)\end{matrix}$it is possible to define the frequency semblance

$\begin{matrix}{{\rho_{a}\left( {f,{p;\overset{\Cap}{\alpha}},\overset{\Cap}{\varphi}} \right)} = {\frac{1}{K\;{\Pi(f)}}{{\sum\limits_{l = 1}^{L}{\mathbb{e}}^{{\{{{- {({\overset{\Cap}{\alpha} + {A{(f)}}})}} + {{\mathbb{i}}{({\overset{\Cap}{\varphi} + {2\;{\pi{({{pf} - {k{(f)}}})}}}})}}}\}}{\delta{({l - l_{0}})}}}}}^{2}}} & (74)\end{matrix}$and write

$\begin{matrix}{{\rho_{a}\left( {{p;\overset{\Cap}{\alpha}},\overset{\Cap}{\varphi}} \right)} = \frac{\int{{{\overset{\sim}{X}\left( {a,f} \right)}}^{2}{\rho_{a}\left( {f,{p;\overset{\Cap}{\alpha}},\overset{\Cap}{\varphi}} \right)}{\mathbb{d}f}}}{\int{{{\overset{\sim}{X}\left( {a,f} \right)}}^{2}{\mathbb{d}f}}}} & (75)\end{matrix}$Further defining

$\begin{matrix}{{\gamma = {{- \frac{\delta}{2}}\left( {\overset{\Cap}{\alpha} + {A(f)}} \right)}}{v = {\frac{\delta}{2}\left( {\overset{\Cap}{\varphi} + {2\;{\pi\left( {{pf} - {k(f)}} \right)}}} \right.}}} & (76)\end{matrix}$yields

$\begin{matrix}{{\rho_{a}\left( {f,{p;\overset{\Cap}{\alpha}},\overset{\Cap}{\varphi}} \right)} = {\frac{1}{K\;{\Pi(f)}}{{\sum\limits_{l = 1}^{L}{\mathbb{e}}^{2{({\gamma + {{\mathbb{i}}\; v}})}{({l - l_{0}})}}}}^{2}}} & (77)\end{matrix}$After some algebra, this simplifies to

$\begin{matrix}{{\rho_{a}\left( {f,{p;\overset{\Cap}{\alpha}},\overset{\Cap}{\varphi}} \right)} = {\frac{\sinh\left( {\overset{\Cap}{\alpha}\;\delta} \right)}{\sinh\left( {L\;\overset{\Cap}{\alpha}\;\delta} \right)}\frac{\sinh\left( {{A(f)}\delta} \right)}{\sinh\left( {{{LA}(f)}\delta} \right)}{\frac{\sinh\left( {L\left( {\gamma + {{\mathbb{i}}\; v}} \right)} \right)}{\sinh\left( {\gamma + {{\mathbb{i}}\; v}} \right)}}^{2}}} & (78)\end{matrix}$considering the hyperbolic sine function with complex arguments. Finallydefining

$\begin{matrix}{{{{\Upsilon(\gamma)} = {\frac{\sinh\left( {\overset{\Cap}{\alpha}\;\delta} \right)}{\sinh\left( {L\;\overset{\Cap}{\alpha}\;\delta} \right)}\frac{\sinh\left( {{A(f)}\delta} \right)}{\sinh\left( {{{LA}(f)}\delta} \right)}\frac{\sinh^{2}\left( {L\;\gamma} \right)}{\sinh^{2}(\gamma)}}}{\Lambda(\gamma)} = {\frac{1}{\sinh^{2}(\gamma)} - \frac{L^{2}}{\sinh^{2}\left( {L\;\gamma} \right)}}},} & (79)\end{matrix}$yields the second order Taylor expansion around v=0 of the frequencysemblance asρ_(a)(f, p

≈Υ(γ){1−Λ(γ)v²}  (80)where γ and v are defined in (76) and focus is directed to the regionaround v=0 for the analysis since that value maximizes the semblance andtherefore, with good estimates of the phase factor, this is the regionof interest. Note, however, that since

≈A(f) in the frequency band of the wavelet coefficients, the attenuationquantity γ is not close to 0 in general and so a similar Taylorexpansion for γ is not made.

Now exploiting the spectral concentration of the wavelet coefficients{tilde over (X)}(a, f) (defined as in (73)) around the center frequencyf_(a), expand the wavenumber, k(f), and attenuation, A(f), around f_(a)using the Taylor series expansion as was done in equation (33). However,for the analysis, retain an extra term in each case, i.e. take theTaylor expansion up to the second order term for k(f) and the linearterm for A(f)

$\begin{matrix}\begin{matrix}{{k(f)} = {{k\left( f_{a} \right)} + {{k^{\prime}\left( f_{a} \right)}\left( {f - f_{a}} \right)} + {\frac{k^{''}\left( f_{a} \right)}{2}\left( {f - f_{a}} \right)^{2}} + {o\left( {{f - f_{a}}}^{2} \right)}}} \\{\approx {{s_{\phi}f_{a}} + {s_{g}\left( {f - f_{a}} \right)} + {\frac{s_{g}^{l}}{2}\left( {f - f_{a}} \right)^{2}}}} \\{{A(f)} = {{A\left( f_{a} \right)} + {{A^{\prime}\left( f_{a} \right)}\left( {f - f_{a}} \right)}}} \\{\approx {A_{0} + {A_{1}\left( {f - f_{a}} \right)}}}\end{matrix} & (81)\end{matrix}$where the quantities s_(φ) etc. above are understood to be the values atf_(a) of each of the corresponding frequency dependent quantities and(·)′ denotes the derivative with respect to f. Using these expansions,it is possible to express v from equation (76) as

$\begin{matrix}{{v \approx {\frac{\delta}{2}\left( {\overset{\Cap}{\varphi} + {2\;{\pi\left( {{pf} - {s_{\phi}f_{a}} - {s_{g}\left( {f - f_{a}} \right)} - {\frac{s_{g}^{i}}{2}\left( {f - f_{a}} \right)^{2}}} \right)}}} \right)}}\mspace{14mu} = {v_{0} + {v_{1}\overset{\Cup}{f}} + {v_{2}{\overset{\Cup}{f}}^{2}}}} & (82)\end{matrix}$where {hacek over (f)}=f−f_(a) is defined as the natural frequencyargument around the band center f_(a) and the coefficients v₀, v₁ and v₂are defined as

$\begin{matrix}{{v_{0} = {\frac{\delta}{2}\left( {\overset{\Cap}{\varphi} + {2\pi\;{f_{a}\left( {p - s_{\phi}} \right)}}} \right)}},{v_{1} = {\delta\;{\pi\left( {p - s_{g}} \right)}}},{v_{2} = {{- \frac{\delta\pi}{2}}s_{g}^{i}}}} & (83)\end{matrix}$It is similarly possible to express γ as

$\begin{matrix}{{\gamma \approx {{- \frac{\delta}{2}}\left( {\overset{\Cap}{\alpha} + A_{0} + {A_{1}\left( {f - f_{a}} \right)}} \right)}}\mspace{14mu} = {\gamma_{0} + \gamma_{ɛ} + {\gamma_{1}\overset{\Cup}{f}}}} & (84)\end{matrix}$with

$\begin{matrix}{{\gamma_{0} = {{- \delta}\; A_{0}}},{\gamma_{ɛ} = {{- \frac{\delta}{2}}\left( {\overset{\Cap}{\alpha} - A_{0}} \right)}},{\gamma_{1} = {{- \frac{\delta}{2}}A_{1}}}} & (85)\end{matrix}$taking into account that γ₀ is the value of γ for the correct estimateof attenuation and γ_(ε) is the perturbation due to the error in theestimate.

With the above approximations, the quantities in equation (79) can beexpressed as

$\begin{matrix}{{{{\Lambda(\gamma)} \approx {\Lambda\left( \gamma_{0} \right)}} = {\frac{1}{\sin\;{h^{2}\left( \gamma_{0} \right)}}\frac{L^{2}}{\sin\;{h^{2}\left( {L\;\gamma_{0}} \right)}}}}{{\Upsilon(\gamma)} \approx {1 - {{\Lambda\left( \gamma_{0} \right)}\left( {\gamma_{ɛ} - {\gamma_{1}\overset{\Cup}{f}}} \right)^{2}}}}} & (86)\end{matrix}$where the terms due to the perturbations, γ_(ε) and γ₁{hacek over (f)},in Λ are not included since it is assumed that these are small and thelatter is already multiplied by similarly small terms in v. Those termsare therefore dropped since it is desired to capture the first ordereffects of the perturbations. Putting this back into the expression (80)for the frequency semblance yieldsρ_(a)(f, p;

)≈1−Λ(γ₀)(γ_(ε)−γ₁{hacek over (f)})²−Λ(γ₀)(v₀+v₁{hacek over(f)}+v₂{hacek over (f)}²)²   (87)The overall semblance can now be computed as in equation (75) by takingthe spectrally weighted average over frequency. However the dependenceof the semblance on the test parameter, p, and the estimated parameters,

and

, is captured to first order by the two complement terms above, thefirst dependent on

(actually since a good estimate of attenuation requires an appropriatevalue of p, there is an implicit dependence on the latter as well), andthe second dependent on p and

. It is therefore possible to write

$\begin{matrix}{{{\rho_{a}\left( {{p;\overset{\Cap}{\alpha}},\overset{\Cap}{\varphi}} \right)} \approx {1 - {\rho_{a,1}^{c}\left( \overset{\Cap}{\alpha} \right)} - {\rho_{a,2}^{c}\left( {p,\overset{\Cap}{\varphi}} \right)}}}{with}} & (88) \\{{{\rho_{a,1}^{c}\left( \overset{\Cap}{\alpha} \right)} = {{\Lambda\left( \gamma_{0} \right)}\frac{\int{{{\overset{\sim}{X}\left( {a,f} \right)}}^{2}\left( {\gamma_{ɛ} - {\gamma_{1}\overset{\Cup}{f}}} \right)^{2}{\mathbb{d}f}}}{\int{{{\overset{\sim}{X}\left( {a,f} \right)}}^{2}{\mathbb{d}f}}}}}{{\rho_{a,2}^{c}\left( {p;\overset{\Cap}{\varphi}} \right)} = {{\Lambda\left( \gamma_{0} \right)}\frac{\int{{{\overset{\sim}{X}\left( {a,f} \right)}}^{2}\left( {v_{0} + {v_{1}\overset{\Cup}{f}} + {v_{2}{\overset{\Cup}{f}}^{2}}} \right)^{2}{\mathbb{d}f}}}{\int{{{\overset{\sim}{X}\left( {a,f} \right)}}^{2}{\mathbb{d}f}}}}}} & (89)\end{matrix}$where {hacek over (f)}=f−f_(a) and {tilde over (X)}(a, f) incorporatesthe attenuation normalization Π(γ) as in equation (73).

In order to complete this computation, the spectral moments are definedas follows:

$\begin{matrix}{{{{f_{\delta} = \frac{\int{{{\overset{\sim}{X}\left( {a,f} \right)}}^{2}\overset{\Cup}{f}{\mathbb{d}f}}}{\int{{{\overset{\sim}{X}\left( {a,f} \right)}}^{2}{\mathbb{d}f}}}}\Delta_{f}^{2}} = {\frac{\int{{{\overset{\sim}{X}\left( {a,f} \right)}}^{2}\left( {\overset{\Cup}{f} - f_{\delta}} \right)^{2}{\mathbb{d}f}}}{\int{{{\overset{\sim}{X}\left( {a,f} \right)}}^{2}{\mathbb{d}f}}}\mspace{34mu} = {\frac{\int{{{\overset{\sim}{X}\left( {a,f} \right)}}^{2}\left( {f - f_{a}} \right)^{2}{\mathbb{d}f}}}{\int{{{\overset{\sim}{X}\left( {a,f} \right)}}^{2}{\mathbb{d}f}}} - \left( f_{\delta} \right)^{2}}}}{{\Gamma_{f}^{3} = {\frac{\int{{{\overset{\sim}{X}\left( {a,f} \right)}}^{2}\left( {\overset{\Cup}{f} - f_{\delta}} \right)^{3}{\mathbb{d}f}}}{\int{{{\overset{\sim}{X}\left( {a,f} \right)}}^{2}{\mathbb{d}f}}}\mspace{31mu} = {\frac{\int{{{\overset{\sim}{X}\left( {a,f} \right)}}^{2}\left( {f - f_{a}} \right)^{3}{\mathbb{d}f}}}{\int{{{\overset{\sim}{X}\left( {a,f} \right)}}^{2}{\mathbb{d}f}}} - {3f_{\delta}\Delta_{f}^{2}} - \left( f_{\delta} \right)^{3}}}},}} & (90)\end{matrix}$representing respectively, the difference between the spectrallyweighted mean frequency and the center frequency f_(a) the spectrumspread (variance) around the spectral mean frequency, and the ‘skew’ ofthe spectrum around the mean frequency. It is then possible to obtainthe semblance complement quantities of equation in terms of the above

$\begin{matrix}{{{\rho_{a,1}^{c}\left( \overset{\Cap}{\alpha} \right)} = {{\Lambda\left( \gamma_{0} \right)}\left\lbrack {\left( {\gamma_{ɛ} - {\gamma_{1}f_{\delta}}} \right)^{2} + {\gamma_{1}\Delta_{f}^{2}}} \right\rbrack}}{{\rho_{a,2}^{c}\left( {p;\overset{\Cap}{\varphi}} \right)} = {{\Lambda\left( \gamma_{0} \right)}\left\{ {\left\lbrack {v_{0} + {v_{1}f_{\delta}} + {v_{2}\left( {f_{\delta}^{2} + \Delta_{f}^{2}} \right)}} \right\rbrack^{2} + {\frac{1}{\Delta_{f}^{2}}\left\lbrack {{v_{1}\Delta_{f}^{2}} + {v_{2}\left( {\Gamma_{f}^{3} + {2\; f_{\delta}\Delta_{f}^{2}}} \right)}} \right\rbrack}^{2} + {v_{2}^{2}(\ldots)}} \right\}}}} & (91)\end{matrix}$where the quantity ( . . . ) is independent of the estimated parametersand so is not described here in any greater detail.

These semblance complements are minimized when the squared quantitiesare zero. Therefore, setting those to zero and substituting thedefinitions of the coefficients from equations (83) and (85) yields thefollowing expressions for the estimates that maximize the semblance:

$\begin{matrix}{{\overset{\Cap}{\alpha} = {{A\left( f_{a} \right)} + {{A^{\prime}\left( f_{a} \right)}f_{\delta}}}}{\hat{p} = {{s_{g}\left( f_{a} \right)} + {s_{g}^{\prime}\left( {f_{\delta} + \frac{\Gamma_{f}^{3}}{2\Delta_{f}^{2}}} \right)}}}{\overset{\Cap}{\varphi} = {2{\pi\delta}\left\{ {{{s_{\phi}\left( f_{a} \right)}f_{a}} - {\hat{p}f_{a}} + {\frac{s_{g}^{\prime}}{2}\left( {\Delta_{f}^{2} - f_{\delta}^{2} - \frac{\Gamma_{f}^{3}f_{\delta}}{\Delta_{f}^{2}}} \right)}} \right\}}}} & (92)\end{matrix}$

Those skilled in the art will recognize that the proposed approach doesnot estimate the parameters

and

by explicitly maximizing the semblance as shown above; rather they werecomputed for each choice of p (and window position t) using equation(53). However, it is possible to assume and even verify numerically thatthe computed values corresponding to {circumflex over (p)} are close tothose satisfying the equations in (92) above. Moreover, it is possibleto refine these computed estimates via local searches to satisfy themaximizing equations above. Therefore, it is possible to use theseexpressions to declare the estimates for the attenuation and the groupand phase slowness. First consider the case where the approximations ofequations (33) and (34) hold good, i.e., when the terms from the nextorder of the Taylor expansion (second order for the wavenumber, firstorder for the attenuation) are relatively negligible over the spectralsupport of the coefficientss_(g) ^(l)Δ_(f) □ s_(g), s_(φ), A′(f_(a))Δ_(f) □ A(f_(a))   (93)with the further understanding that

$f_{\delta},{\frac{\Gamma_{f}^{3}}{\Delta_{f}^{2}} < \Delta_{f} < {f_{0}.}}$Substituting these approximations leads to the estimates

$\begin{matrix}{\left. {\hat{A}\left( f_{a} \right)}\leftarrow{\overset{\Cap}{\alpha} \approx A_{0}} \right.\left. {\hat{s}}_{g}\leftarrow{\hat{p} \approx s_{g}} \right.{\left. {\hat{s}}_{\phi}\leftarrow{\hat{p} + \frac{\overset{\Cap}{\varphi}}{2\pi\; f_{a}\delta}} \right. = {{{s_{\phi}\left( f_{a} \right)} + {\frac{s_{g}^{\prime}}{2f_{a}}\left( {\Delta_{f}^{2} - f_{\delta}^{2} - \frac{\Gamma_{f}^{3}f_{\delta}}{\Delta_{f}^{2}}} \right)}} \approx {s_{\phi}\left( f_{a} \right)}}}} & (94)\end{matrix}$

Note that these approximations hold good even when the conditions in(93) are not satisfied so long as the odd moments are negligiblecompared to the bandwidth, i.e.,

$f_{\delta},{\frac{\Gamma_{f}^{3}}{\Delta_{f}^{2}}{{\bullet\Delta}_{f}.}}$That holds when the spectrum represented by {tilde over (X)}(a, f) iswell centered on f_(a). This is the case when the spectrum of the modeof interest is relatively flat over the effective support of theanalyzing wavelet at scale a.

Bias Correction

When neither of the conditions stated above hold good, i.e., when thederivatives of the group slowness and the attenuation are notinsignificant and when there is significant spectral asymmetry, then theestimates from equation (94) are biased. The resulting bias can behandled in at least two ways.

Frequency correction is motivated by the bias in both the attenuationand the group slowness estimates in (92) containing a linear term thatis the linear Taylor expansion of the corresponding quantity evaluatedat f_(a)+f_(δ). Assuming that this dominates the other terms, the biascan be handled by declaring the estimates in (94) to actually refer tothe frequency f_(a)+f_(δ) instead of the center frequency, f_(a). Thephase slowness estimate, which shows smaller bias, can be appropriatelycorrected by estimating it as if the center frequency was f_(a)+f_(δ).This correction is logical since this is the effective center frequencyfor the frequency band of the processed coefficients. Note that f_(δ)can be readily calculated from the data. While simpler to implement,this can change the frequency support of the computed dispersion curveestimates, often shrinking it relative to the processing frequency band.

Slowness correction is another approach which explicitly estimates thebias terms in (92) and corrects for them. This process includesestimating the spectral moments f_(δ), Δ_(f) ² and Γ_(f) as well as thederivatives of group slowness and attenuation. The former can beestimated from the data. The latter can be obtained once the groupslowness and attenuation dispersion curves have been estimated. Theseare the biased dispersion curves and are liable to be noisy. However, byusing smoothing operators to handle the noise and exploiting thecondition that the bias in the derivatives is relatively small comparedto the bias of the estimates themselves (due to smoothness of dispersioncurves), a first order bias correction can be achieved. It is possibleto further refine this estimate by another iteration if necessary.

FIG. 10 shows an example of bias in the group slowness estimation in theregion of large slowness variation and spectral asymmetry and itscorrection by both methods on a synthetic data example. The grey lineshows the true value. Note that the support of the dispersion curveafter frequency correction shrinks to 8.3 kHz under frequencycorrection. Also note that the phase estimates are relatively unaffectedfor this case.

Dispersion Curve Extraction

The techniques described above in connection with attenuation andslowness (group and phase) analysis at one scale of the array CWT forone mode can be extended to other scales (and corresponding centralfrequencies) to obtain an estimate of slowness for this mode as afunction of frequency, i.e., the dispersion. However, such an approachcan be computationally costly. An alternative approach exploits thedependency across scales arising from the continuity and smoothness ofphysical dispersion curves, the similar continuity across scales of thetime frequency localization of the mode and the relationship of phaseand group slowness. In particular, given the phase and group slowness ata certain scale, corresponding to ω_(i) where ω_(i) could refer eitherto the nominal center frequency, 2πf_(a) or to the spectral meanfrequency 2π(f_(a)+f_(δ)), it is possible to deduce the followingrelationship between the phase correction and the phase and groupslowness at the frequency, ω_(i+1), using a Taylor series expansion

$\begin{matrix}{{\phi\left( \omega_{i + 1} \right)} = {{\omega_{i}{s_{\phi}\left( \omega_{i} \right)}} + {\frac{\omega_{i + 1} + \omega_{i}}{2}{s_{g}\left( \omega_{i} \right)}} + {o\left( {\omega_{i + 1} - \omega_{i}} \right)}}} & (95)\end{matrix}$The above relationship can be used to limit the search over the phaseand use the continuity of the modes in time and slowness across scalesto likewise limit the search over the time and slowness domains. Thisnot only results in the code being more efficient, but by employing thiscontinuity to track across scales it is also possible to obtain improvedaccuracy and robustness in the estimates. The described embodimenttherefore utilizes the approach of initializing the estimation byevaluating the time frequency map at a reference sensor (usually thelast) where the best separation of the modes might be expected. On thismap, a starting value of the frequency (scale) is chosen correspondingto maximum energy (dominant mode). A full search of the EPRT (over thecandidate moveouts p) is carried out for this scale using time windowlocations around the time of maximum modulus. This yields estimates forthe phase and group slowness at this frequency. The technique then“marches up” the frequency using the estimates from the previous stepsto constrain the search in the present step. The marching (tracking) isan iterative increase (or decrease) which may be terminated based ontesting the semblance against a threshold. After marching up, thetechnique “marches down” from the initial starting point, similarlyterminating by comparing the semblance to a threshold (the order andtiming of marching up versus down is not critical, and no limitation isimplied). Note that it is now possible to extract dispersion curves formore than one mode using the same approach if they do not overlap in thetime-frequency plane. After tracking the strongest mode, the process isre-initialized with the next strongest peak corresponding to the secondmode and the process repeated to extract the corresponding dispersioncurve. This process can be iterated until the significant modes that aredisjoint on the time-frequency plane have been tracked and thecorresponding dispersion curves extracted. In this regard, it is helpfulto use the data association algorithm described above to track andcorrectly label the modulus peaks. FIG. 7 shows an example of using dataassociation to track multiple modes on the time frequency plane.

Another use of data association could be to label the modes acrosssensors, fit a line to the estimated time locations as a function ofsensor and use the slope there from to initialize the search over p overa small neighborhood of the fitted slope. This could lead to furthersavings in computation. However there may be an issue of whether thedata association is always robust, e.g., when there are interferingmodes nearby. This modification effectively uses the relatively lightmodulus method described previously to constrain the search whilerefining the estimates obtained therefrom. This process is illustratedin FIG. 4 and FIG. 5. Note that FIG. 5 also illustrates collecting thecoefficients at a certain frequency (scale) into an array indexed bytime and sensor.

An example of application of the EPRT method to real oilfield data isshown in FIG. 11, where it is compared to the existing matrix pencilmethod. It will be appreciated that the new method yields more stableand better estimates for both slowness and attenuation.

The above-described approach works for cases where the modes do notoverlap in the time-frequency plane, though they could overlap in timeor frequency. In practice, this should cover a majority of cases ofinterest. However, in some instance there is significant overlap ofdistinct modes in the time frequency plane. Another variation topossibly address that case is to use the reassigned scalogram describedabove to obtain a modified time-frequency (time-scale) map showingbetter definition and separation of the modes. It is then possible toapply the method described here, provided the reassignment operationremoves the overlap. Methods to address the case of overlapping modesthat are not separated in the reassigned scalogram will be describedbelow.

Workflow

The EPRT workflow used in extracting the slowness and attenuationdispersion curve estimates is shown in FIG. 12. The steps are asfollows:

1. Start with the waveforms collected by an array of at least twosensors.

2. Compute the continuous wavelet transform (CWT) coefficients for eachwaveform, generating a time-frequency (time-scale) map for a desired setof frequencies or scales (as shown in FIG. 4).

3. Select a reference sensor (preferably one which shows greatestseparation between the modes of interest); on the corresponding CWT map,identify and track the separated modes across frequency

a) Identify and obtain the time locus of each mode using a feature onthe CWT map. Examples of such features include the peak of the absolutevalue or modulus (possibly after data fitting), onset defined as whenthe modulus first exceeds a threshold relative to the peak, locus ofphase satisfying a given value near the modulus peak etc. FIG. 4 showsthe peak locations.

b) Use data association to track these time loci across frequencies andthereby label the individual modes across frequency (as shown in FIG.7).

Optionally repeat above for the waveform at each sensor. Then do dataassociation across sensors of the modes labeled across frequency toidentify and label them across sensors (as shown in FIG. 5).

1) For each of these labeled modes:

-   -   a) Initiate the dispersion curve extraction at the frequency        corresponding to the greatest energy for that mode.    -   b) Collect the CWT coefficients at this initial frequency for        the sensors into an array of data indexed by time and sensor (as        shown in FIG. 5).    -   c) Pick time windows (the window length is decided based on the        frequency) centered on a set of time positions, t, in a local        neighborhood of the time locus of the mode corresponding to the        initial frequency.    -   d) Pick a set of test moveouts, p, corresponding to the expected        range of group slowness. If the labeling has been done across        sensors in the optional step above, this set can be constrained        by fitting a straight line through the corresponding time loci        of the mode of interest for the sensors, computing its slope and        choosing the set of test moveouts to be in the local        neighborhood of this computed slope.    -   e) For each of the test moveouts, shift and align the array of        coefficient data corresponding to the moveout.    -   f) For each of these aligned arrays, compute estimates of the        exponential parameters, phase        and attenuation,        , as described previously, for the data in each of time windows        applied to the shifted data.    -   g) Use these estimated exponential parameters to compute the        EPRT semblance according to equation (66). Build up the map        ρ(t, p) and locate its maximum peak. The corresponding moveout        is declared as the estimate {circumflex over (p)}.    -   h) Use the estimate {circumflex over (p)} as well as the        corresponding values of        and        to compute the estimates of group and phase slowness and        attenuation at this frequency as per equation (94).    -   i) Alternatively, generate a map from the numerator of the        semblance, corresponding to the energy of the projection, locate        its peak and use the corresponding values of p,        and        to generate our dispersion estimates.    -   j) The dispersion curve estimation is then continued by marching        up the frequency axis; moving onto the next higher selected        frequency, collecting the CWT coefficients corresponding to it        in a data array and repeating steps (4a) through (4i) above. The        modification is that the search is now constrained over test        moveouts, p, to lie in a certain neighborhood of the estimated        value of {circumflex over (p)} at the previous frequency value.        The phase estimates are also constrained to lie within a        neighborhood of a predicted value based on previous estimates as        per equation (95). If the computed phase estimate lies outside        this interval, set the estimate to equal the interval endpoint        closest to it.    -   k) This marching up process is terminated when the highest        available frequency is reached, or when the maximum computed        semblance falls below a specified threshold.    -   l) The process is repeated by marching down in frequency        starting with the one below the initial frequency.    -   m) When done, a dispersion curve for this mode results. As a        final step, apply the bias correction if desired. This could be        either of the following:        -   i) Frequency correction: Estimate the energy spectrum            weighted frequency mean for the CWT coefficients in the            window corresponding to the maximizing t and p values at            each of the selected frequencies. Assign the group slowness            and attenuation to these mean frequencies and re-compute the            phase slowness estimates using these mean frequencies            instead of the center frequencies as described before. This            correction could also be performed in the marching steps            above.        -   ii) Slowness correction: Compute the energy spectrum            weighted spectral moments—frequency mean, frequency variance            (bandwidth squared) and skew (third moment) as per            equation (90) for the CWT coefficients in the window            corresponding to the maximizing {circumflex over (t)} and            {circumflex over (p)} values at each of the selected            frequencies. In addition, calculate the local slope after            smoothing of the uncorrected group and attenuation estimates            as a function of frequency. Use these derivative estimates            along with the spectral moments to apply corrections to the            slowness and attenuation estimates as per equation (92).    -   2) When done, the dispersion curve estimates are obtained over a        subset of the processing frequency band for the phase and group        slowness as well as the attenuation for each of the identified        and labeled modes. These modes are assumed herein to be        sufficiently separated on the time frequency map of the received        to permit this labeling.

While the invention is described through the above exemplaryembodiments, it will be understood by those of ordinary skill in the artthat modification to and variation of the illustrated embodiments may bemade without departing from the inventive concepts herein disclosed.Moreover, while the preferred embodiments are described in connectionwith various illustrative structures, one skilled in the art willrecognize that the system may be embodied using a variety of specificstructures. Accordingly, the invention should not be viewed as limitedexcept by the scope and spirit of the appended claims.

1. A method of providing a dispersion curve for acoustic data comprisingthe steps of: using an analyzer unit to receive a time series ofacoustic data associated with a plurality of sensors having differenttransmitter-receiver spacings; generate a plurality of time-frequencyrepresentations from the received time series; identify a characteristicfeature in the time-frequency representations; select at least onefrequency; define a kernel comprising: estimating a time location of thecharacteristic feature at the selected frequency; associating acorresponding characteristic feature for ones of the spacings at theselected frequency; fitting at least a segment of a line through atleast some of the time locations of the associated characteristicfeatures, plotted against the corresponding spacings; estimating groupslowness at the selected frequency based at least in-part on anindication of fitted line segment slope; estimating a phase differencebetween sensors; estimating an attenuation between sensors; using theestimated phase difference and group slowness to estimate a phaseslowness; and repeat the kernel for other frequencies of interest,obtain a dispersion curve comprising the group slowness, the phaseslowness, and the attenuation as a function of frequency, and providethe dispersion curve in tangible form.
 2. The method of claim 1including the further step of using the analyzer unit to generate atime-frequency map by applying at least one step selected from the groupincluding: a wavelet transform to the time series for each sensor; aWigner Wille transform to the time series at each of the sensors; and ashort time Fourier transform to the time series at each of the sensors.3. The method of claim 1 wherein said characteristic feature is at leastone of: a contour of peaks with respect to time of a complex modulus ofa time freguency map; and a contour of a complex modulus of a timefrequency map that equals a given threshold.
 4. The method of claim 1including the further step of using the analyzer unit to apply areassigned scalogram operation to a time-frequency map.
 5. A method ofproviding a dispersion curve for acoustic data comprising the steps of:using an analyzer unit to: receive a time series associated with aplurality of sensors; generate a time-frequency representation for eachsensor; select a reference sensor; identify a characteristic feature inthe time frequency representation for the reference sensor; estimatetime locations across frequency of the characteristic feature; select aninitial frequency at which energy is greater than a threshold value overthe time locations; define a kernel comprising: collecting timefrequency coefficients at the selected frequency for the sensors into arepresentation of data indexed by time and sensor; selecting timewindows centered on a set of positions in a local neighborhood of thetime location at the selected frequency; selecting a set of testmoveouts for ones of the time windows corresponding to expected range ofgroup slowness; shifting and aligning the representation of coefficientdata in the window corresponding to the test moveout and windowpositions; computing estimates of phase difference and attenuationacross the representation of shifted and aligned data for the testmoveouts and window positions; using the estimated phase difference andattenuation to perform a modified stacking operation on therepresentation of shifted and aligned data for the test moveouts andwindow positions to produce an output representation for a plurality ofvalues of time position and moveout; finding a maximum peak of theoutput representation and using the corresponding moveout to estimatethe group slowness; using a corresponding value of the computed estimateof the phase difference along with the group slowness to generate anestimate of the phase slowness; setting a corresponding value of thecomputed estimate of attenuation as an attenuation estimate; iterativelyincreasing frequency from the selected initial frequency comprisingselecting a higher frequency; constraining the moveout based on theprevious computed value of the group slowness; repeating the kernel forthe selected higher frequency until a final highest frequency is reachedor when the maximum computed semblance falls below a specifiedthreshold; iteratively decreasing frequency from the selected initialfrequency comprising selecting a lower frequency; constraining themoveout based on the previous computed value of the group slowness;repeat the kernel for the selected lower frequency until a final lowestfrequency is reached or when the maximum computed semblance falls belowa specified threshold; obtain a dispersion curve comprising the groupslowness, the phase slowness, and the attenuation as a function offrequency; and provide the dispersion curve in tangible form.
 6. Themethod of claim 5 including the further steps of using the analyzer unitto: estimate a bias correction; and the bias correction to thedispersion curve.
 7. The method of claim 6 wherein the bias is at leastone of: a frequency bias; and a slowness and attenuation bias.
 8. Themethod of claim 5 including the further step of using the analyzer unitto generate a time-frequency map by applying a step selected from thegroup including: a wavelet transform to the time series at each of thesensors; a Wigner Wille transform to the time series at the of saidsensors; and a short time Fourier transform to the time series at eachof the sensors.
 9. The method of claim 5 wherein said characteristicfeature is at least one of: a contour of peaks with respect to time of acomplex modulus of a time frequency map; and a contour of a complexmodulus of a time frequency map that equals a given threshold.
 10. Themethod of claim 5 wherein the window length is decided based on thefrequency.
 11. The method of claim 5 wherein the output representationis an output map indicative of at least one of: energy; and semblance.12. The method of claim 5 further including the step of using theanalyzer unit to apply a reassigned scalogram operation to atime-frequency map.
 13. Apparatus for dispersion extraction for acousticdata comprising: a sonic logging tool configured to generate a timeseries of acoustic data associated with a plurality of sensors havingdifferent transmitter-receiver spacings; an analyzer unit configured to:generate a plurality of time-frequency representations from the timeseries; identify a characteristic feature in the time-frequencyrepresentations; select at least one frequency; define a kernel,including: estimating a time location of the characteristic feature atthe selected frequency; associating the corresponding characteristicfeature for ones of the spacings at the selected frequency; fitting atleast a segment of a line through at least some of the time locations ofthe associated characteristic features; estimating group slowness at theselected frequency based at least in-part on slope of the fitted linesegment; estimating a phase difference between sensors; estimatingattenuation between the sensors; using the estimated phase differenceand group slowness to estimate the phase slowness; and repeating thekernel for other frequencies of interest, obtaining a dispersion curvecomprising the group slowness, the phase slowness, and the attenuationas a function of frequency, and storing the dispersion curve in memory.14. The apparatus of claim 13 wherein the analyzer unit is furtheroperable to generate a time-frequency map by applying at least one of: awavelet transform to the time series for each sensor; a Wigner Willetransform to the time series at each of the sensors; and a short timeFourier transform to the time series at each of the sensors.
 15. Theapparatus of claim 13 wherein said characteristic feature is at leastone of: a contour of peaks with respect to time of a complex modulus ofa time frequency map; and a contour of a complex modulus of a timefrequency map that equals a given threshold.
 16. The apparatus of claim13 wherein the analyzer unit is further operable to apply a reassignedscalogram operation to a time-frequency map.
 17. Apparatus fordispersion extraction for acoustic data comprising: a sonic logging toolconfigured to generate a time series associated with a plurality ofsensors; an analyzer unit configured to: generate a time-frequencyrepresentation for ones of the sensors; select a reference sensor;identify a characteristic feature in the time frequency representationfor the reference sensor; determine time locations across frequency ofthe characteristic feature; select an initial frequency at which energyis maximum over the time locations; define a kernel, including:collecting time frequency coefficients at the selected frequency forones of the sensors into a representation of data indexed by time andsensor; selecting time windows centered on a set of positions in a localneighborhood of the time location at the selected frequency; selecting aset of test moveouts for ones of the time windows corresponding toexpected range of group slowness; shifting and aligning therepresentation of coefficient data in the window corresponding to onesof the test moveout and window positions; computing estimates of phasedifference and attenuation across the representation of shifted andaligned data for ones of the test moveouts and window positions; usingthe estimated phase difference and attenuation to perform a modifiedstacking operation on the representation of shifted and aligned data forones of the test moveouts and window positions to produce an outputrepresentation for values of time position and moveout; finding amaximum peak of the output representation and using a correspondingmoveout to estimate group slowness; using a corresponding value of thecomputed estimate of the phase along with group slowness to generate anestimate of phase slowness; setting a corresponding value of thecomputed estimate of attenuation as an attenuation estimate; iterativelyincreasing frequency from the selected initial frequency comprisingselecting a higher frequency; constraining the moveout based on theprevious computed value of the group slowness; repeating the kernel forthe selected higher frequency until a final highest frequency is reachedor when the maximum computed semblance falls below a specifiedthreshold; iteratively decreasing frequency from the selected initialfrequency comprising selecting a lower frequency; constraining themoveout based on the previous computed value of the group slowness;repeating the kernel for the selected lower frequency until a finallowest frequency is reached or when the maximum computed semblance fallsbelow a specified threshold; obtaining a dispersion curve comprising thegroup slowness, the phase slowness, and the attenuation as a function offrequency; and storing the dispersion curve in memory.
 18. The apparatusof claim 17 wherein the analyzer unit is further operable to: determinea bias correction; and apply the bias correction to the dispersioncurve.
 19. The apparatus of claim 18 wherein the bias is at least oneof: a frequency bias; and a slowness and attenuation bias.
 20. Theapparatus of claim 17 wherein the analyzer unit is further operable togenerate a time-frequency map by applying at least one of: a wavelettransform to the time series at each of the sensors; a Wigner Willetransform to the time series at each of the sensors; and a short timeFourier transform to the time series at each of the sensors.
 21. Theapparatus of claim 17 wherein said characteristic feature is at leastone of: a peak with respect to time of a complex modulus of a timefrequency map; and a contour of a complex modulus of a time frequencymap that equals a given threshold.
 22. The apparatus of claim 17 whereinthe window length is decided based on the frequency.
 23. The apparatusof claim 17 wherein the output representation is an output mapindicative of at least one of: energy; and semblance.
 24. The apparatusof claim 17 wherein the analyzer unit is further operable to apply areassigned scalogram operation to a time-frequency map.